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Investigation of the Temperature Distribution 
and Thermal Stresses in a Hypersonic Wing 
Structure 


MARTIN A. GOLDBERG? 
Bell Aircraft Corporation 


SUMMARY 


The transient temperature distribution and thermal stresses 
in an idealized stiffened skin stringer panel isolated out of a multi- 
cell wing are calculated for a typical hypersonic aircraft. The 
heat flux into the system, resulting from aerodynamic heating, is 
considered to vary sinusoidally with time and to remain constant 
with position, being equal to the heat flux entering the section 
at a station midway between the stiffeners. Conduction through 
the structural members is analyzed, while conduction, convec 
tion, and radiation through the air trapped within the structure 
is neglected 

Temperature and thermal stress distributions are presented 
in dimensionless parameters for a typical flight plan 

The simplifying assumptions employed in the analysis are dis- 
cussed and their effect on the temperature and stress distribution 
is indicated 

The thermal stresses present in a section are affected most 
seriously by the ratio of web depth to the span of the skin. As 
this ratio increases, the thermal stress in the web passes through 
a maximum, then decreases, whereas the stress in the flange in- 
creases continuously. Eventually the maximum thermal stress 
in the section occurs in the flange. The variation of thermal 
stress caused by changing the other geometric ratios does not 
appear to be as significant as that mentioned above. Curves 
illustrating their effects are presented in terms of dimensionless 
parameters 
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SYMBOLS 


modulus of elasticity 
maximum flux 
dimensionless parameter, 27,/) 
w/O, 

dimensionless parameter, 2b,,./b 
dimensionless parameter, 7,./27 

length of skin between webs 

half the length of the web 

specific heat of material, B.t.u./lb.°F 

x~W/M/2 

(1-—%*x*)V¥ M/2 

heat-transfer coefficient, B.t.u./ft.? hr.°F 

thermal conductivity of material, B.t.u. in./ft.* hr.°F 
Laplace transform notation 

flux, B.t.u./ft.? hr 

dimensionless flux, |(b./2 
temperature at station 7 
equilibrium temperature 
initial temperature, °F 
temperature ratio at station 7, 


RV M/2 


dimensionless position coordinate 


dimensionless position coordinate 


(1 —R)VM/2 

thermal diffusivity, k/pc, 

coefficient of thermal expansion 

time 

exposure time 

dimensionless time for analog solution, a6 


dimensionless time for Laplace solution, 


dimensionless constant, ).£)/2kt, 
density 

thermal stress 

dimensionless thermal stress, o;/a@Et, 
thickness of skin 

thickness of web 


merical computations; and to Miss E. Hautala for her help in the Subscript I relates to the flange 


Preparation of the numerical computations Subscript II relates to the web 








982 JOURNAL OF THE 





od 











TEMPERATURE OF THE SKIN 


CHORDWISE STATIONS 








<C| | Loo 


TYPICAL SUPERSONIC WING STRUCTURE 
CONSIDERED IN ANALYSIS 


Chordwise temperature distribution for a supersonic 
wing at small angles of attack 





























Fic. 1 


INTRODUCTION 


_ STRUCTURAL DESIGN of high-speed aircraft is 
strongly influenced by consideration of the effects 
of aerodynamic heating. Theory predicts boundary- 
layer temperatures of appreciable magnitude for 
supersonic vehicles, and these predictions are borne out 
with a reasonable degree of accuracy by experiments. 

A problem of considerable interest is the determi- 
nation of the temperature distribution and thermal 
stresses in multiweb structures. Various solutions of 
this problem have been proposed in the past, based on 
different simplifying assumptions. It is the purpose 
of the present paper to examine the available solutions, 
analyzing what is inherently implied in the assump- 
tions made, and to present a new solution based on 
heat input assumptions that are believed to correspond 
more closely to the actual physical situation. 

Prior to this investigation, some experimental work 
on built-up structural sections had been done by 
Barzelay and Boison' with the data on temperature 
and thermal stress distributions tabulated in graphical 
form. 
tered in actual supersonic flight were not simulated; 


However, conditions similar to those encoun- 


therefore, the usefulness of the data is limited. 

Torda and Hoff" have treated a problem in which a 
However, in their work 
they solved the problem in two parts—in the first part 
the heat loss from the cap to the web was assumed 
negligible and in the second the temperature rise of 
the web at the junction point of the cap and web was 
an exponential function of time, as determined by the 
temperature rise of the cap at the same point from 
the first part of the problem. The analysis used by the 
authors yielded good results for the configuration con- 
sidered since the cap was very heavy and the web very 
thin, the ratio of cap area to web area being 20 to 1. 
However, for a large number of practical cases of in- 


“T”’ section was analyzed. 
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terest, the geometry of the section would conflict with 
the assumptions, and the resulting solution could leaq 
to erroneous results. As a result of neglecting the 
conduction from the cap into the web the cap tem 
perature will appear to rise more rapidly than in the 
actual case. Consequently, the web, since it is sub 
jected to a higher heating rate, will also exhibit a tem 
perature rise somewhat greater than in the true physical 
situation. The net result is the prediction of lower 
thermal stresses than actually occur. Therefore, the 
solution yields optimistic values. Torda and Hoff also 
considered a constant heat-transfer coefficient and 
boundary-layer temperature. According to the avail 
able theory, the heat-transfer coefficient depends upon 
the speed, altitude, and temperature of the skin; there. 
fore, the variation of the coefficient with time may be 
quite pronounced. The boundary-layer temperature 
is known to vary approximately as the square of the 
Mach Number, and the possible variation is consider- 
able. The consequences of using the above assumption 
are to force the flux entering the skin to be a maximum 
immediately upon launching the craft, and to decrease 
as time increases. In the true case, the flux is quite 
small at the time of launching, building up to a maxi 
mum as the velocity increases, after which time it be- 
gins to decrease either because the equilibrium tem- 
perature is reached or the craft decelerates. 

In the light of what has been said it can be seen that, 
with the method of solution presented by Torda and 
Hoff, the thermal stresses are predicted to occur earlier 
in the flight than is actually the case. Also, the struc- 
ture is actually heated more gradually than is indicated 
in the solution presented by Torda and Hoff, which 
will probably result in lower values of peak thermal 
stresses than are estimated by the authors. 

Three years after the above paper was published, 
Pohle and Oliver!” presented a more exact treatment of 
the problem. By using the same numerical values of 
Torda and Hoff, they bore out the latters’ results and 
conclusions within 1 per cent accuracy. In this paper, 
as in the aforementioned one, only constant heat- 
transfer coefficient and constant boundary-layer tem- 
perature were treated. Here, as before, the work is 
based on heat transfer from a boundary layer of con- 
stant temperature, and may lead to overdesign as 4 
result of the flux input into the system. However, 
the general method presented is extremely useful within 
the limitations of the assumptions, and is employed in 
the present analysis. 

This paper will treat a case similar to that described 
by Pohle and Oliver, but with different heat input 
assumptions. In an effort to improve upon the as 
sumptions in the latters’ analysis, a method was formu- 
lated whereby the effects of varying the heat-transfer 
coefficient and boundary-layer temperature with time 
could be taken into consideration. The flux entering 
the structure is calculated on the basis of a predeter- 
mined flight plan, hence permitting a close association 
of this analysis and the performance of the craft. 
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TEMPERATURE 


ANALYTICAL TREATMENT OF THE PROBLEM 


A supersonic wing flying at a small angle of attack 
chordwise distribution, as 


Considering a symmetrical airfoil of 


will have a temperature 
shown in Fig. 1. 
uniform material, and at a small angle of attack, the 
heat transferred to the upper surface will be the same 
is that transferred to the lower surface. In this case, 
the mid-plane of the wing can be considered to be per- 
fectly insulated. In order to simplify the analysis 
somewhat further, the station under investigation is 
chosen at the mid-span and mid-chord. Since tem- 
perature does not vary appreciably in the spanwise 
and chordwise direction in this region, the problem can 
be treated by considering unidirectional heat flow. 
\sa result of the station chosen, any effects due to the 
curvature of the surfaces can be neglected. Using the 
ibove assumptions, the section being examined can be 
isolated. The idealized section is shown in Fig. 2. 

The thicknesses of the materials utilized in aircraft 
construction are generally small compared to their 
other dimensions. Therefore, the temperature §gra- 
dients through the thicknesses are neglected; only those 
along the length of the members are considered. 

The problem is simplified still further by assuming 
that the physical and thermal properties of the material 
Only the aerody- 
through the 
convection, and 


are independent of temperature. 
namic heating effects and conduction 
member are treated; 
radiation through the air trapped within the structure 


The errors introduced by these assump- 


conduction, 


are neglected. 
tions are examined in another section. 

Typical of many aircraft design calculations, the 
computation of the aerodynamic heating effects re- 
quires a knowledge of the expected performance or 
fight plan. The the boundary-layer 
temperature versus time can be calculated from the 
fight plan, using the procedure established in reference 
13. With the aid of the heat balance equation from 
the same reference, a transient skin temperature cal- 
culation can be performed for a point midway between 
the stiffeners yielding a relationship of skin temper- 
Since the inner surface 


variation of 


ature as a function of time. 
is perfectly insulated, the flux entering the skin can be 
estimated by equating it to the heat absorbed by the 


body. The net result of all the preliminary calcula- 
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Idealized skin stringer section 
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Fic.4. Simplified representation of member under consideration 


tions is to obtain a curve of the heat input with respect 
to time. A typical curve of this nature is shown in 
Fig. 3. 
at which cooling begins, and the slope of the curve 


The maximum magnitude of the flux, the time 


depend upon the exact flight plan. 
In order to utilize the flux input in the analysis, the 
actual flux is approximated by the relationship 


g = Fy sin (1, 6¢)6 (1 


where £; is the maximum flux, 6, is the duration of flux 
input that is, from time zero to cooling—and @ 1s 
time. 

As can be seen from Fig. 2, the section is symmetrical 
about two axes; therefore, only one quarter need be 
considered and can be represented as shown in Fig. 4. 

The general differential equations governing the flow 
of heat in the section shown in Fig. 4 is derived from 
the basic heat balance in an element. For any station 
in the skin, where 6, < x < b,/2 + 6, and y = 17;, the 


unidimensional heat conduction equation is 

(O°t;/Ox?) + (g/Rkrs) = (1, @) (Ot1, 08) (2 

In the web, where 0 < x < b,, vy = (7,/2), the conduc 
tion equation is 

O7t1;/Ox? = (1) a@) (Of, 06 (.3) 


Introducing the dimensionless variables 


A aad v a tt 
iia . x= “= 
[(b,/2) + 6b, ]?’ (6,/2) + 5 
[(b,/2) + b,;PE: . b,, T 
= . R = ’ S = 
. tok, ' (b,/2) + b, Qt 


and substituting these in Eqs. (2) and (3), we obtain 
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O71 _ WO Ou; f cs (d2uy*/dx?) + [gM /(p? + M?)] = puy* uy, (9 
+ g sin = for R<*X<1 (4) Ig f sia : 
Ox? 6, fal 
d7uy1* dx? = purn* = oe 10) 
Oni : Ou11 F <xe< - s 
av2—sog lor OS4S8h (9) where 
The associated boundary conditions for Eqs. (4) and M = 2/€ 
(5) are = ws 
lhe transformed boundary conditions are 
| i i x= 1 Ou;/OX = 0 
x=R uy = Un dl, duy*/dX¥ = 0 | 
(1/.S) (Ou;/Ox) = Ony,/Ox (6 x=R uy* = un" 
gc (0 Ou1/0x = 0 ’) (1/S) (duy*/dX) = duy* ™ 
a= 6 uy = 1 =O duyy*/dX¥ = 0 
ur = | ; ? ' ve 
Letting r> = p, the general solutions for Eqs. (9) and 
If wu = u(x; 6), then the Laplace transform of uv will (10) are 


be denoted by u*, where 


: I gM 
. : uy* = Ccosh rx + D sinh rx + [ _ 
u* = u*(X; p) = | e ” u(x; 9) d9 (7) p (p? + M?) | 
vi) 12 
The Inversion Theorem‘ then states that = 
fork <. 2% < I 
I eae pe 
u(x; 9) = is { e” u*(x; p) dp (S) uyy* = A cosh rx + B sinh rX¥ + (1/p 13 
ont J 4-i ; 


for0 < x < R, where A, B, C, and D are arbitrary con 


over a suitable contour in the complex p-plane. 
stants of integration. 


The Laplace transform of Eqs. (4) and (5) yields the 
ordinary differential Eqs. (9) and (10) for the trans- Solving for the constants of integration using the 
forms of uw; and uy. boundary conditions given in Eqs. (11), they are 


7 gM tanh r(1 — R) 
7 p(p? + AM?) cosh rR [tanh r(1 — R) + S tanh rR] 
B=0 15 


C —gMS tanh rR cosh r ) 

= = = = (10 
p(p? + M?*) cosh r(1 — R) [tanh r(1 — R) + S tanh rR] 

gMS tanh rR sinh r 


D= “ ~ = 
p(p? + M*) cosh r(1 — R) [tanh r(1 — R) + S tanh rR} 


Substituting the values for the arbitrary constants in Eqs. (12) and (13), they become 
a | [ gM | gMS tanh rR cosh r(1 — £) 
= = bs é . (18 
(p? + M?) p(p? + M*) cosh r(1 — R) [tanh r(1 — R) + S tanh rR] 


l 4 gM tanh r(1 — R) cosh rf 
p p(p? + M?) cosh rR [tanh r(1 — R) + S tanh 7R] 





un* = 





In order to obtain the equations for u; and “1, the Inversion Theorem is applied. Therefore, with the use of Eq 


(8), Eqs. (18) and (19) become 


g GMS rte tanh uR cosh u(1 — 2) e*? dy m 
(wy — 1) = (1 — cos M@) — : [ : é = = : = 2U 
M 2mi J,-i2 WwW? + M?) cosh u(1 — R) [tanh w(1 — R) + S tanh uR] 
1) GM prt tanh u(1 — R) cosh yu e”? dy (91 
uit — = = = = <1 
se 2mi J,-i2 Ww? + M*) cosh uR [tanh w(1 — R) + S tanh pR] 


where r = yu, p = y, and yp = Vy. If we take the same branch of the radical in the numerator and denominator 0! 
the integrands, it is seen that the numerator and denominator are single-valued analytic functions. Therefore, the 
appropriate contour to use is the familiar one given in Fig. 29 of reference 4. If the integral over the contour ap 
proaches zero as the radius tends to infinity, then the evaluation of Eqs. (20) and (21) is equivalent to the evaluation 
of the residues of the integrand. The proof that the integral over the curved portion of the contour tends to zere 
follows the procedure generally used with problems of this nature and is given in reference 4. 
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The integrand of Eq. (20) has poles at 


y = 0 


y = +1M 

y= Be ae s fal with m = 0, 2, 4 se 
1(1 — R)? — 

y = —¢, where ¢, is given by tan ¢,(1 — R) + Stan ¢,R = 0 


while the poles of Eq. (21) are 


y =0 
y = +1M | 
(m + 1)*x° ; (23 
y=- = with m = 0,2,4.... 
tR? 
y = —?¢,” where ¢, is given by tan ¢,(1 — R) + Stan oR = 0 


In the limit, the line integrals given by Eqs. (20) and (21) are equal to 277 times the sum of the residues at the 
poles given by Eqs. (22) and (23), respectively. The proof of this is given in reference 4 and need not be given here 
After the necessary manipulations and simplifications have been performed, the temperature ratios are 
RS 
M(1 — R + RS) 
2S tan o,R cos ¢,(1 — Loe 
,(,? + M2) cos ¢,(1 — R) [C1 — R) sec? ¢,(1 — R) + RS sec? ¢,R] 
2S tan VR cos Y(1 — ®)e~** 
V(¥++ iM) (1 — R) sin VU — R) [tan VL — R) + Stan YR] 
(C; cosh g cos g + C2 sinh g sin g) cos ad on 


— 4 (1 — cos V@) — au | T 


. V/ 


dD, 
vl (1 — R) 2 tan ¢,(1 — R) cos ¢,%e 
Wy — 1 = g. = “re oe = - = - — “pret "hes 
Wi -—- R+ RS) o,(o,1 + AM?) cos ¢,R[(1 — R) sec? o,(1 — R) + RS sec? ¢,R] 
r 5 zy Z% . . . ; . 
2 tan Z(1 — R) cos Zl e (C3; cosh d cos d + C; sinh d sin d) cos M@ | - 
_ = = — . (25 
Z(Z* + M?)R sin ZR [tan Z(1 — R) + Stan ZR] Ds, | 
where 
(n + 1)x 
jhe p> BP (26 
n=0, 2,4 2(1 — R) 
. al (n + 1l)r 
Z= = (97 
re 9 = 
n=O, 2, 4 2R 
(sinh 27 sinh z cos z + sin 2v sin z cosh 2) 2.S(sinh? v + sin* v) cosh z cos 2 
( ; = = : . me + . ? : : : (2S 
S(cos? v + sinh? v) S(cos? v + sinh? v 
: (sinh 2v sin z cosh z — sin 2v sinh z cos 2) 2.S(sinh? v + sin’ v) sinh zg sin 2 
_= % + " P . (29 
S(cos? v + sinh? v) S(cos? v + sinh? v) 
: 2(sinh? z + sin? z) cosh v cos v S(sinh 22 sinh v cos v + sin 22 sin v cosh v 
C,; = r - 5 (30 
; 2(sinh? z + sin? z) sinh v sin v S(sinh 2s sin v cosh v — sin 22 sinh v cos v 
CG, = : sw : (31) 
a 2 
(ae habs (cos? s + sinh? s) (sinh? v + sin? v) | sinh 22 sinh 2v + sin 22 sin 27] _,., 
D, = 2M — (sinh? z + sin* 2) + : oe a - rer (32) 
Pe (cos? v + sinh? v 2S cos? v + sinh* v 
' 2(sinh? v + cos” v) (cosh? z — cos* 2) S(sinh 2v sinh 22 + sin 2v sin 22) = 
D, = 2M*| (cosh? v — cos? v) sone “ + ee (+309) 
2S?(sinh? z + cos? 2) 2.S?(sinh? z + cos? 2 
Having obtained the temperature distribution can be calculated from the equations’ 
through the member, we have finally to evaluate the 6, = Bi — un (34 
thermal stresses. The thermal stresses in a thin plate o,= Bi — uy (35 








YS6 
where 
] | l" a 
By = — > uy dt + | 1x 
als es tat eed Fal 


(36) 


The use of the above equations for the determination 
of the thermal stresses is contingent on the assumptions 
that the temperature gradient through the member is 
symmetric with respect to the centroid of the section 
and that the edges, x = 1, are free. Since the skin is 
relatively weak in bending out of its plane of maximum 
moment of inertia, it is reasonable to consider that the 
edges, £ = 1, are free for the thermal stress calculation. 


NUMERICAL RESULTS 


The analysis presented has been applied to seven 
different configurations to determine the effects of vary- 
ing the geometrical ratios with the same flux input. 
The calculations were performed on an electronic dif- 
ferential analyzer, being based on a nondimensional 
flux equal to 1.25 sin (70/0.075). This flux is the 
amount necessary to cause the skin of a structural “T”’ 
section, having the dimensions 27,/b, = 0.07, (2b,,/b,) = 
0.5, and 7,/27, = 0.4, to reach its equilibrium tem- 
perature midway between stiffeners at the nondimen- 
sional time, 6 = 0.075. A typical temperature dis- 
tribution obtained in this manner is given in Fig. 5 for 
three points in the skin and four in the web. The 
geometrical proportions employed to obtain the curves 
are indicated on the Figures. 

In addition to the above calculations, the temper- 
ature distribution was determined for a section with 
half of the total heat input that was used in the previous 
calculations. The data are given in Fig. 6. 

Fig. 7 shows the temperature gradients through 
the ‘‘T’’ section for various times. This is a cross plot 
taken from a temperature distribution curve and is 
utilized to determine the thermal stresses. 


DISCUSSION 


Discussion of Results 


Comparing the temperature distribution curves ob- 
tained, the effects of varying the geometry can be ascer- 
tained. In order to facilitate this comparison, curves 
are presented of the maximum temperature difference as 
a function of time for the four parameters investigated. 

The effect of increasing the web depth can be deter- 
mined from Fig. 8. Keeping three parameters constant 
and varying the ratio Xk, the maximum temperature 
difference increases as 6, is increased. Examining 
Fig. 8, it can be seen that, when & is equal to 1.5, the 
maximum possible temperature difference is rapidly 
being approached. For all practical purposes, any 
increase in R beyond the value of 1.5 will have a negli- 
gible effect on the temperature difference. 

Since the thermal stresses are based on the difference 
between the weighted average temperature and the 
temperature at the point in question, the effect of in- 
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creasing K on the maximum thermal stresses is cop 
siderably different than on the temperature differences 
Fig. 9 illustrates the effect of increasing 6, on the mayj 
mum thermal stress. As noted before, increasing p 
increases the temperature difference at a decreasing 
rate, whereas it has the opposite effect on the weighted 
average temperature. The existence of these two 
opposite effects yields a peak thermal stress in the web 
at a particular value of XK. As XK is increased beyond 
this point, the stress in the web decreases and th 
stress in the flange increases. Eventually, the mayj 
mum thermal stress will be in the flange. 

The variation in temperature difference caused by 
increasing the thickness of the web, all other parameters 
being constant, is quite small. This can be seen from 
Fig. 10. The values of S employed in the analysis wer 
chosen to represent the practical range of 7,,. 27,. 

Fig. 11 shows the effect of varying the web thickness 
on the maximum thermal stress in the web. The stress 
decreases as 7, increases, and since the temperature 
difference also decreases in a like manner, eventually a 
point will be reached at which the thermal stresses in 
the flange are greater than those in the web. This is not 
significant, however, since it appears to be out of the 
range of practicability. 

The ratio of the skin thickness to the semispan of 
the skin, 77, can be seen from Fig. 12 to govern the 
overall temperature level. As it was just shown that 
the ratio S had little effect on the temperature differ- 
ence, increasing 7, serves to decrease the temperature 
level of operation. This is understandable, as in- 
creasing 7, merely increases the heat capacity of the 
skin. Therefore, a longer period of time is necessary for 
the skin to reach equilibrium temperature. 

Although the change in temperature level of opera- 
tion is quite pronounced with an increase in 7,, the 
effect on the thermal stresses is considerably less, as 
seen in Fig. 13. Actually, Fig. 13 is based on a con 
stant value of S; however, increasing /7 decreases S 
since K is constant. If this variation had been taken 
into consideration, the variation in the maximum 
thermal stress would be less pronounced. As an ap- 
proximation, consider changing S from 0.5 to 04, 
which results in a 3.5 per cent increase of the maximum 
thermal stress, whereas changing /7 from 0.07 to 0.09 
results in a decrease of the maximum thermal stress of 
8.0 per cent. In order to change H from 0.07 to 0.09 
with R = 0.5, the corresponding change in S would be 
approximately from 0.5 to 0.4. Hence, it can be seen 
that the effects are in opposite directions yielding 2 
less pronounced variation in the maximum. thermal 
stress with // than indicated. 

The difference in temperature level of operation due 
to the variation of exposure time is quite pronounced, 
and can be seen readily from Fig. 14. Any flux dur- 
ation less than that necessary to cause the skin to be at 
equilibrium temperature results in higher temperatures 
throughout the member at any instant. The maximum 
temperature differences are also higher; however, the 
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maximum thermal stress in the web is considerably 
lower. 

The variation of thermal stress with time is shown in 
Fig. 16. It is interesting to note the change in stress 
that occurs at the junction, x = 0.324, from tension 
to compression. Although the stress in the flange is 
small compared to that in the web, it is of sufficient 
magnitude to warrant investigation of the stability of 
the flange, since it is usually critical in compression 
under normal conditions. As the ratio R increases, 
this consideration will gain in importance as mentioned 


previously. 


Sources of Error 


The assumptions used in this analysis induce a cer 
tain amount of error in the determination of the tem- 
perature gradients and thermal stresses. It is the pur- 
pose of this section to examine the assumptions and any 
possible sources of error that might exist and to show 
their effects on the calculations. The assumed condi- 
tions will be examined in the order of their estimated 


importance. 


Conduction Through the Air 


Conduction through the air trapped within the air- 
foil, which was neglected in the analysis, can have a 
considerable effect upon the temperature gradient in 
the member. The conductivity of air, while increasing 
with increase of temperature, is practically constant 
with pressure, until the mean free path of the molecules 
approaches the size of the container, at which time it 
varies directly with pressure. Therefore, as long as a 
craft is within the atmosphere, the air contained 
within the structure will transmit heat. The effect of 
this additional heat being carried from the skin to the 
stringer is to help decrease the temperature gradient in 
the web and to reduce the peak thermal stresses. The 
neglect of this added source of heat input to the stiffener 
tends to cause the solution to yield conservative results. 


Effect of Neglecting Radiation 


In the analysis presented, radiation between the 
flanges and the webs were neglected. The familiar 
equation governing this mode of heat transfer is given by 

g = oF ,F.(T\4 — T>') (37) 
Recognizing that heat is radiating to the web from 
four individual sources at the same temperature, the 
flux obtained through the use of Eq. (87) is one-fourth 
of the total flux entering the web. An estimate of the 
view factor, Fy, to employ in the above equation can be 
obtained from Fig. 4-10 of reference 10. 

The heat transferred by radiation can, of course, be 
decreased by coating the interior of the structure with 
a low emissivity surface, such as gold. However, it is 
obvious that it would be most advantageous to assist 
this mode of heat transfer. Therefore, coating the 
inside surfaces with a high emissivity substance, such 
as carbon black, would facilitate the heat transfer by 
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radiation between the skin and the webs. The overalj 
effect of increasing the radiant heat flow would be ty 
decrease the temperature gradients in the web, thereby, 
reducing the thermal stresses in the section. 

From the above discussion, it can be seen that 
neglecting radiation yields larger temperature gradients 
hence, the predicted thermal 


than actually exist; 


stresses are conservative. 


Flux Input Constant with Position 


In the analysis presented, the variation of the fluy 
entering the system with position was neglected. The 
actual heat entering the section at some time early ip 
the flight will have a distribution as shown in Fig. 17 
The flux was assumed to be constant with a magnitude 
It can 
be seen from the diagram that a quantity of heat less 
than that which exists in the true physical situation js 
Consequently, the 


equal to that shown at the edges of the skin. 


permitted to enter the structure. 
temperature rise in the web will be less rapid than js 
the actual case, resulting in thermal stresses of a higher 
order of magnitude. As time progresses, the difference 
between the assumed and true values of the flux will 
diminish. The overall effect of considering this type 
of flux distribution is to yield a conservative estimate of 
the thermal stresses present in the section. 


Junction Points 


At the junction of the ideal stiffener and skin perfect 
contact was assumed. With the exception of integrally 
stiffened skins, this does not exist. 
be a small surface film or air space around the area 


There will always 
where the two members are joined. The presence of 
the film has a decided insulating value; consequently, 
the heat being conducted from the skin to the web will 
be somewhat retarded, yielding a larger temperature 
drop between the skin and the web. If the members 
are welded, the presence of the foreign material will 
also act to delay the heat from entering the stiffener t 
some degree. The net results of these types of inherent 
insulations are to raise the temperature in the skin and 
decrease the temperature in the web. Consequently 
the thermal stresses will occur a little earlier in time 
than predicted by the theory; the peak values will be 
substantially larger. It is hoped that this error will be 
offset by the conservative errors introduced in the 


previous assumptions. 


Effect of Convection Within the Airfoil 

Free convection of the air trapped within the airfoil 
was neglected on the basis of being small compared t 
the other modes of heat transfer in the system. Unfor 
tunately, at the present time there is not sufficient infor 
mation available to estimate the heat transferred b) 
free convection in the section under consideration 


Large Angles of Attack, or Unsymmetrical Construction 


For the case when the airfoil is not constructed sy™- 
metrically, either in gages or type of construction, or !s 
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flying at large angles of attack, the simplification lead- 
ing to the analysis of one-quarter of the section cannot 
be applied. In that case, half the section must be 
analyzed taking into consideration the flux entering 
both the top and bottom surfaces of the section. 


Leading Edge 


The analysis presented in the paper is based on uni- 
directional heat flow. Therefore, the solution obtained 
can only be employed in regions remote from the lead- 
ing edge. In the area adjacent to the leading edge, a 
two-dimensional analysis must be employed. 


Number of Divisions Used in the Electronic Differential 

Analyzer Solution 

The accuracy of the results obtained when using the 
electronic differential analyzer is dependent upon the 
number of linear differential equations utilized to re- 
place the original partial differential equation. As 
in the case of all numerical solutions, the finer the net- 
work used, the more accurate will be the results. 

The structural ‘‘T’’ section analyzed in this paper 
was divided into fifteen sections. The skin was divided 
into eleven elements in order to ensure that the junc- 
tion of the skin and web would occur at a half element, 
thereby simplifying the equations employed. 

The error introduced by using the number of sec- 
tions employed in the analysis is of the same order of 
magnitude as that encountered when calculating the 
heat flux entering the section. This error is within 
engineering accuracy. 


CONCLUSIONS 

(1) The overall temperature level of a section is 
governed by the ratio of the skin thickness to the 
semispan of the skin. 

(2) The maximum temperature difference is greatly 
dependent upon the ratio of the web depth to the span 
of the skin, increasing at a decreasing rate as the ratio 
is increased. 

(3) The ratio of half the web thickness to the skin 
thickness does not affect the temperature distribution 
in the section to any appreciable extent for the ratios 
investigated. 

(4) A shorter duration of flux input than is required to 
reach equilibrium temperature will cause every point 
in the member to be at a higher temperature at any 
particular time. However, since the total flux input 
is less, the member does not reach the same temperature 
level. 

(5) The thermal stresses present in a section are 
affected most seriously by the ratio of web depth to the 


span of the skin. For small values of this ratio the 
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maximum thermal stress occurs in the center of th 
web, whereas at larger values the peak thermal streg 
occurs in the skin. At some intermediate value th 
maximum thermal stress has a peak with reference t 
this ratio. 

(6) Increasing the thickness ratio, 7,/27,, decreases 
the maximum thermal stresses in the web while increas 
ing those in the skin. However, this variation js ; 
too pronounced. 

(7) Increasing the ratio of the skin thickness to sem; 
span of skin decreases the thermal stresses in the wet 
and increases the maximum thermal stress in the flange 
Small variations in this parameter do not cause sign; 
ficant changes in the peak thermal stresses. 

(8) Any variation in the duration of flux input, rela 
tive to that necessary for the skin to reach equilibrium 
temperature, results in a considerable decrease of th 
maximum thermal stress present in a section. 


REFERENCES 


1 Barzelay, M. E., and Boison, J. C., Investigation of Stre 
Due to Thermal Gradients in Typical Aircraft Structures, NACA 
RM 51K06, January, 1952. 

2 Bosworth, R. C. L., Heat Transfer Phenomena, John Wik 
and Sons, Inc., New York, 1952 

> Brull, M. A., Thermal Stresses in Stiffened Structures 
craft Corporation, Unpublished Report, 1955. 

*Carslaw, H. S., and Jaeger, J. C., Conduction of Heat 
Solids, Oxford University Press, London, 1948 

5 Churchill, R. V., Introduction to Complex 
Applications, McGraw-Hill Book Company, Inc., New York 
1948 

6 Churchill, R. V., Modern Operational Mathematics in Eng 
neering, McGraw-Hill Book Company, Inc., New York, 1944 

7 Dusinberre, G. M., Analysis of Heat Flow 
McGraw-Hill Book Company, Inc., New York, 1949 

8 Jakob, M., Heat Transfer, Vol. I, John Wiley and Sons 
Inc., New York, 1949. 

’ Korn, G. A., and Korn, T. M., Electronic Analog Computer 
McGraw-Hill Book Company, Inc., New York, 1952 

1 McAdams, W. A., Heat Transmission, McGraw-Hill Book 
Company, Inc., New York, 1954. 

11 Peirce, B. O., A Short Table of Integrals, Ginn and Company 
New York, 1920. 

12 Pohle, F. V., and Oliver, H., Temperature Distribution and 
Thermal Stresses in a Model of a Supersonic Wing, Journal of 


8-16, Januar) 


Variables and 


Numerical 


the Aeronautical Sciences, Vol. 21, No. 1, pp 
1954. 

13 Johnson, H. A., Rubesin, M. W., Sauer, F. M., Slack, E.G 
and Possner, L., A Design Manual for Determining the Therma 
Characteristics of High Speed Aircraft, Army Air Force’s Technica 
Report No. 5632, September, 1947. 

14 Timoshenko, §S., and Goodier, J. N., Theory of Elastict' 
McGraw-Hill Book Company, Inc., New York, 1951 

16 Hoff, N. J., Structural Problems of Future Aircraft, Proceed- 
ings, Third Anglo-American Aeronautical Conference at Brighton 
England, 1951, The Royal Aeronautical Society, pp. 77-11 
Appendix I: ‘Aerodynamic Heating and Thermal Stresses, 
by Paul Torda and N. J. Hoff, pp. 103-110. 








simp 
recen 
pote! 
the 1 
an in 
lar t 
assur 
In 
ary-li 
solve 
simp] 
tion 1 
the le 
lamin 
with 
regiol 
copte 
chord 
lamin 
so th 
to th 
by co 
Th 
soluti 
lamin 
have 
Rec 
' es 


ter of the 


mal stress 


value th, 
ference t 


decreases 


le increas 
ion is} 


SS to sem 


n the web 


the flange 


LUSE Signi 


iput, rela- 
juilibrium 
ise of the 


1 of Stre 
res, NACA 
John Wile 
es, Bell A 
Heat 


of 


riables and 


New York 
cs in Eng 
‘k, 1944 
Teat Flow 
and Sons 
Computers, 
-Hill Book 
Company 
bution and 
Journal of 
, January 
ack, E. G., 
1e Therma 
; Technica 
Elastictty 
}, Proceed- 
Brighton 


», 77-114 


Stresses, 


I 
t 





Some Examples of Laminar Boundary-Laver 
Flow on Rotating Blades’ 


NICHOLAS ROTT? ann WILLIAM E. SMITH? 


Cornell University 


SUMMARY 


Laminar flow over a cylindrical blade rotating steadily in an 
incompressible fluid has been studied by Fogarty. The results 
* his investigation of this three-dimensional boundary-layer- 

w problem are here applied to a number of specific examples 
involving both positive and negative pressure gradients in the 
hordwise direction 

In these examples, the chordwise component of flow in the 
yundary layer is the solution of the two-dimensional problem 
treated by Falkner and Skan 


letermined by a linear partial differential equation involving the 


The spanwise-flow component is 


chordwise component This equation is separated into two 


rdinary differential equations by means of appropriate simi- 


larity transformations 





The results of a numerical integration indicate how the span- 


vise component is affected by the chordwise pressure gradient 


INTRODUCTION 


} AMINAR BOUNDARY-LAYER PROBLEMS concerning 


L rotating blades have become soluble, after certain 
simplifications, on the basis of several investigations 
recently published. First, it is necessary to know the 
potential flow about rotating blades. Sears! has given 
the notably simple solution for the potential flow on 
an infinite cylinder rotating about an axis perpendicu- 
lar to its generators. Subsequently, a blade will be 
assumed to be cylindrical. 

In the case of a rotating blade, the laminar bound- 
ary-layer problem has been attacked by Fogarty’ and 
solved in several special cases. The most important 
simplifying assumption made by Fogarty (in addi- 
tion to the usual boundary-layer assumptions) is that 
the length of chord on which the development of the 
laminar boundary layer occurs is small in comparison 
In 


region including most of the span of a propeller, heli- 


with its distance from the center of rotation. a 
copter rotor, or possibly even a turbine blade, the 
chord length of 
laminar boundary-layer separation is sufficiently small 


between leading edge and _ point 


so that this assumption is reasonable. Confinement 
to this small-chord region makes the problem soluble 
by comparatively simple means. 

The purpose of the present paper is to give further 
Solutions to Fogarty’s problem. In two-dimensional 
laminar boundary-layer theory, Falkner and Skan® 
have given a set of exact solutions whose fundamental 
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property is that the velocity profiles are similar at any 
chordwise station. It is suggested, then, that corre- 
sponding solutions are obtainable on rotating blades. 
Actually, for the ‘spanwise’? boundary-layer flow in 
this three-dimensional case, similarity is generally im- 
possible to obtain, as is shown later; however, a solution 
is found by combination of two profiles, each similar, 
with variable ‘‘weight.”’ 

For the purposes of this investigation some discus- 
of the 
are needed and are given first. 
tions show the effect of rotation on laminar boundary 


sions and extensions aforementioned papers 


The subsequent solu- 


layers for several cases of accelerated and decelerated 
flow. 
FLow 


(1) THe POTENTIAI 


Consider an infinite cylinder of arbitrary cross sec 
tion rotating with constant angular velocity, w, about 
an axis normal to its generators. Introduce a rotating 
coordinate system, X, Y, Z (see Fig. 1), the 
Z-axis is the axis of rotation and the cylinder generators 


where 
are parallel to the }-direction. The velocity compo- 
nents L’, |’, IW of an incompressible fluid are given in the 
rotating system by the following solution due to Sears.! 
Let gil X, Z) be the two-dimensional potential of the 
flow past the same cylinder placed in a parallel stream 


of unit speed in the positive Y-direction. Then 


U = w Very V = WY] 2x WW wl giz (1.1 


The solution is readily verified since (a) the diver- 
gence of the velocity vector vanishes, (b) its curl vector 
has only one nonvanishing component which is parallel 


to the Z-axis, viz., 


(c) the boundary conditions on the cylinder are satis- 
fied, and (d) the conditions at infinity are fulfilled. 
Far from the cylinder ¢:(X, Z) becomes the potential 
of a parallel stream of unit velocity —.e., 


gi(X, Z) ~ X + a constant 1.2) 
In this case, Eqs. (1.1) become 
U = wY Vv = —w(X + a constant) W=0 
(1.3) 


which are really the velocity components in the rotating 


system when the fluid is absolutely at rest. Further- 
more, if the Z-axis is to be the axis of rotation, the con- 


stant in the expression for I’, Eq. (1.3), must be zero. 
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Cartesian and boundary-layer coordinates in the rotating 
system. 


Fic. 1 


This eliminates the indeterminacy in |’ of Eq. (1.1), 
for ¢: has to be normalized so that, far from the cylin- 
der, the constant in Eq. (1.2) vanishes. 

The velocity components of the ‘‘cross flow,” 
W, are given by Eq. (1.1) exactly in the form which is 
also expected by means of a rather naive two-dimen- 
sional analysis. They increase proportionally with 
Y, whereas the flow V is found to be 
independent of the spanwise coordinate. Thus, for 
large Y, V < U and V < W ultimately. A parallel 
shift of the axis of rotation affects only |’, and this by 


U and 


‘ 


‘spanwise”’ 


an additive constant. 

In subsequent applications to the laminar boundary- 
layer solutions, the velocity distribution of the Falkner 
and Skan 
sidered. The chordwise velocity component outside 
the boundary layer, ™, is taken proportional to the 
coordinate x (see Fig. 1), measured along the surface, 


two-dimensional potential flow’ is con- 


raised to a constant power m—1.e., 


Uy = wY(x/c)™ (1.4) 


where c is a characteristic chord length for which “, = 
wY. The value of the potential g, along the surface is 


gi = [c/(m + 1)]) &/o™!' + C (1.5) 


In two-dimensional flow a velocity distribution of 
the type of Eq. (1.4) may be realized by an incom- 
pressible flow past an infinite wedge of half-angle 
am/(m + 1), but is found also in the vicinity of a 
wedge-shaped leading edge of any closed body. The 
latter interpretation may be applied to rotating blades. 
If the common origin of coordinates X and x occurs at 
the leading edge, and if the wedge shape is symmetrical 
about the X-axis, then 


X = x cos mm/(m + 1) (1.6) 


Moreover, for properly shaped bodies, velocity dis- 
tributions of this type can occur at a sharp leading 
edge, even though symmetry about the X-axis is 
lacking. Thus, for generality, Eq. (1.6) may be re- 
placed by 
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X = x cos 7 17 


with arbitrary y. In cases of negative m, a ‘‘negative 
wedge angle appears to be a rather artificial interpre. 
tation; but discussion of this case of retarded flow is of 
special interest. Therefore, it is often advantageoys 
to consider the velocity distribution of Eq. (1.4) to be 
produced by some artificial means over a flat portion 
of a surface inclined at an angle y to the X-direction. 
The spanwise component v; outside the boundary 


layer, according to Eqs. (1.1) and (1.5), is 


v1, = w{[c/(m + 1)] (v/c)™t+! — 2x cosy + C} (18 


The constant C depends upon the axis of rotation. Ip 
particular, there is always an axis of rotation for which 
v, vanishes at the leading edge—i.e., for which C = () 
In the case of a biconvex profile, for example, it can be 
shown that the chordwise position of the proper axis 
is (1 — m) half-chords from the mid-chord toward 
the leading edge. Such problems have been investi- 
gated by Graham.* The value of C for a given axis of 
rotation is the only quantity (necessary for the bound- 
ary-layer calculations) which cannot be deduced from 
a given (x) alone; the shape of the whole body and 
(in principle) the whole potential field g:(X, Z) must 


be known. 


(2) THE BOUNDARY-LAYER EQUATIONS 


In addition to the usual boundary-layer assumptions, 
Fogarty” has introduced two simplifications. The 
first restriction, already mentioned, is to chord lengths 
which are small compared with the distance from the 
axis of rotation; in other words, one considers asymp- 
totic solutions for large radii. The second assumption 
is that the normal to the body surface is everywhere al- 
most parallel to the axis of rotation. 

For the present investigation, it is proposed to inter- 
pret the flow (at least tentatively) as that of a stream 
Hence, it seems desirable to abandon 
It is found that it is pos- 


past a wedge. 
Fogarty’s second restriction. 
sible to generalize Fogarty’s equations for an arbitrary 
angle y between the surface normal and the axis of 
rotation, if only small variations of the angle y are 
The latter restriction is in- 

It may be recalled that 


admitted over the surface. 
significant for wedge flow. 
small curvature of the surface and thus small rate of 
change of y is also a necessary assumption of the usual 
boundary-layer theory. For rotating blades, additional 
“higher-order” terms due to curvature arise, expressing 
combined effects of curvature and rotation. The 
equations have been discussed in full generality by 
Mager.’ 

The derivation of the boundary-layer equations for 
y nonvanishing and nearly constant is not presented 
here, as it follows exactly Fogarty’s procedure. The 
results are given in boundary-layer coordinates on the 
rotating surface (see Fig. 1): x, y, 2 are measured 
chordwise and spanwise along the surface and normal 
to it, respectively; u, v, w are the corresponding velocity 
components in the rotating system; 1; and % are the 
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values of « and v outside the boundary layer. The 
equations are” 
UUy + WU, = Uy, + VU (2.1) 
uv, + WU, + 2w COS yu = Uy, + W-: (2.2) 
Ur +w, = 0 (2.3) 
with the boundary condition 
u=v=w=0 for 3 =O (2.4) 
“=m, Ver, forz = of 
The only difference from Fogarty’s equations is the 


appearance of the factor cos y. 

As a consequence of Fogarty’s first restriction (to 
small chords), the Eqs. (2.1) and (2.3) for the chord- 
wise flow are exactly the same as for a two-dimensional 
boundary layer. No effect of rotation is found in the 
chordwise flow, which is reduced to a known case. 
Subsequently, if « and w are known, the spanwise 
component, v7, can be determined from Eq. (2.2), which 
isa linear inhomogeneous differential equation in 2. 

At chordwise sections located at different spanwise 
positions, the chordwise boundary-layer developments 
are similar, with “,; being proportional to wy. It fol- 
lows from the general similarity properties of laminar 
boundary layers that the thickness, 6, of the layer 


2) whereas the shearing stress 7, ~ 


varies as (wy) 
wy 

The spanwise flow component outside the layer, 7, 
has been found to be independent of y in the previous 
section. Eq. (2.2) for v inside the layer can be made 
independent of y by putting 
u' = u/wy, v’ = v/we, 


2.5) 


¢ = 2/6, 6 = V vC/wy 


72 


By use of the continuity Eq. (2.3), write Eq. (2.2) in 


the form 


*. 
uv, — a | uy dz + 2w cos yu = UU, + Ww, (2.6) 
After introduction of the variables in Eq. (2.5), one 
obtains 
uy’? + v,,' 


os 
ft / ’ Se ‘ 
We — Ue | uz’ dg + 2 cos yu’ = 


which is independent of y and w. This means that the 
magnitude of the spanwise component of the boundary- 
layer profile is also independent of y but the thickness 
of the spanwise profile varies as (wy)~‘!?)—i.e., in the 
The ratio 
‘u varies, therefore, as 1/y along the span. Con- 
sequently, it is not advantageous to discuss the effects 


ol rotation on the boundary layer by considering the 


same manner as for the chordwise profile. 


nondimensional 
IS tO use v/wx or v As 
’ Wis of the order of x/y, examination of vy/ux may be 


fatio vu. The proper way to make v 


wc, which is independent of y. 


desirable. 

A property of the solution of Eq. (2.2) may lead to 
Substantial simplifications, especially for the examples 
Which will be treated subsequently. The linear in- 
homogeneous equation for v, Eq. (2.2), has two inhomo- 


F 
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geneous terms. It is possible to split v into two parts 


Vv =U. + UY (2.8) 


which fulfill, respectively, the following differential 
equations: 
(2.9) 


= MN, T Wo, 


(2.10) 


UUyy + WY, = —2w COS YU + WWp,, 


The boundary conditions to be fulfilled by zy, and 2, are, 


forz = 0, 
% = % = 0 
and forz = ~, we have 
Vas = Uns = Yaz = 1 = Q (2.11) 
These properties of v, and v, at infinity cause Eqs. (2.9) 
and (2.10) to yield 
U1 (Uar)1 = Ul, (2.12) 
U1(Up,)1 = —2w COS yy (2.13) 
According to Eq. (1.1), #1, = w¢i,. Therefore, from 
Eq. (2.12), 

Yq, = w¢r (2.14) 
is the proper boundary condition for v, at 2 = 
Likewise, from Eq. (2.13), it follows that for z = 

Up, = —2w cos yx (2.15) 


Thus, at the limit of the boundary layer, uy, and % as- 
sume the values of the two terms which appear in Eq. 
(1.1) for V. 

For the 
equation for u, to be treated subsequently, the corre- 


“similar’’ solutions of the boundary-layer 


sponding solutions for vy, and v, each have the similarity 
property. This property, which makes a solution pos- 
sible, is not possessed by the total solution for v. 

The sole effect of a parallel shift of the axis of rota- 
tion on the boundary conditions, Eqs. (2.4) is a change 
of 2% by a constant. Consider two solutions of Eq. 
(2.2) with the boundary conditions v 
v1 + wC. The difference between the two solutions, 
v., has the boundary condition at z = 


= vy, and v,. = 


v.. = (2.16) 


ve, = wl = const. 
and evidently satisfies the homogeneous part of Eq. 


(2.2)—+.e., 


to 
~I 


UVer FF Weg = Wee: (: 
The complete solution after a shift can be written 


(2.18) 


with the differential Eqs. (2.9), (2.10), and (2.17) and 
the boundary conditions (2.14), (2.15), and (2.16) to 
be fulfilled. 

Eq. (2.17) is identical to the ‘‘spanwise’’ boundary- 
layer equation for an infinite cylinder at yaw, as given 
Prandtl® and Sears,’ and has the same boundary 


by 
Originally, in the present case, ™ varies 


conditions. 
proportional to wy, but by use of the variables defined 
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in Eq. (2.5), this variation is taken care of and the case 
is reduced to that of the yawed infinite cylinders. 
Thus, the spanwise profiles are identical in both cases; 
only 6 behaves differently. The thickness is constant 
along the infinite cylinder at yaw, while it varies pro- 
portionally to y~''”) on a rotating blade. The change 
in v due to the shift of the axis of rotation has been 
reduced to a known case 

The spanwise boundary layer on a yawed cylinder 
has been determined by Cooke* for the Falkner and 
Skan velocity distribution, Eq. (1.4). The profiles 
are known and have been found by Rott and Crabtree’ 
to be fairly independent of m if reduced to the proper 
scale. Therefore, with no loss of generality, the sub- 
sequent investigation is restricted to one value of C in 
Eq. (1.8), viz., C = 0. 


(3) THE SOLUTION FOR FALKNER AND SKAN VELOCITY 
DISTRIBUTIONS 
The boundary-layer equations, (2.1) to (2.4), are 
to be solved for the velocity distributions given by 
Eqs. (1.4) and (1.8). Determination of the chordwise 
boundary layer leads to the two-dimensional problem 
solved by Falkner and Skan.* 
Let 
u/u, = f'(¢) 
where 
¢ = 32/A = 2/5F(x/c) = V [(m + l)wy]/20c X 


/ 
(x/c)\™@-DI2 g = sV [(m + 1)u,]/2xv (3.2) 


From the continuity equation, (2.5), the normal ve- 


locity component is 


/ 
w= —V [((m + l)myv]/2x {f + 
[((m — 1)/(m + 1) ]¢f"} 


Eq. (2.1) produces an ordinary differential equation 


for f: 
ff +f" + Bl — f’”) =0 (3.4) 
where B = 2m/(m + 1) (3.5) 
The boundary conditions are 
P seu a for¢é = 0 } cdr 
o~ ; de (3.6) 
y =] forg¢ = of 


The functions /’(¢), which give the shapes of the ve- 
locity profiles, have been computed and tabulated for 
different values of m by Hartree." 

Determination of the spanwise boundary layer de- 
pends upon the form of the expression for the spanwise 
velocity outside the boundary layer, Eq. (1.8), with 
C = 0. 
solution of Eq. (2.2) for the spanwise boundary-layer 


Following the discussion in Section (2), the 


component is split into two parts, v, and %, fulfilling 
Eqs. (2.9) and (2.10). 
Assuming “‘similar’’ solutions for both parts—i.e., 


profiles depending upon ¢, one may write 


VY = [we/(m + 1)] (w/c)™+'g(¢) (3.7) 
g 
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_ >) 


% = —2wcos yx h(f) (3.8 


From Eqs. (2.9) and (2.10), two ordinary differentia] 
equations for g and hare obtained: 


+i —- a= -2 (3.9 
h” + fh’ — (2-—6)f'h = -—(2 — B)f’ 5.10 


Both equations are linear and inhomogeneous; / and /’ 


are known functions of ¢ from Hartree’s solution 
Eq. (3.10) becomes homogeneous by putting h — 1 =} 
k” + fk’ — (2 — B)f'k = 0 3.1] 


The boundary conditions are 


=h=0, or k= 1 fors = 0 } 


5 > 19 
« ) 
g= h = 1, or R= 0 iors = { 


The special case of the flat plate, 8 = 0, has been 


treated by Fogarty.? In this case, the homogeneous 
parts of Eqs. (3.9) and (3.10) are identical, so that 
they can be joined into one equation. Generally, the 
homogeneous parts differ in the coefficient of the third 
term, and the functions g and /: (or k) have to be de- 
termined by solutions of different equations. The total 
solution, v = v% + %, is obtained by superimposing the 
profiles g and / with “weights” variable along the 
chord, as they are multiplied by different powers of x 
Eqs. (3.7) and (3.8). 

, 


Eqs. and (3.11) 
Eq. (3.12), have been solved by step-by-step calculations 


(3.9) with boundary conditions, 
It is advantageous to divide the range of the inde- 
intervals which are smaller 


pendent variable ¢ into 
in the neighborhood of ¢ = 0 where the functions g and 
k vary more greatly, and to use larger intervals for 
large values of ¢. This is accomplished by changing 
the independent variable to 

a=In(1+ ¢) (3.13 
and by using equal increments of a. The functions 
g and h: are plotted in Fig. 2 for different values of 3 
Table 1 shows the values g’(0) and h’(0) which are 
necessary to express the spanwise component of the 
shearing stress. 

Another decomposition of the velocity components 
is employed by Fogarty and is appropriate here. Con- 
sider components “* and v* which are parallel and 
perpendicular, respectively, to the resultant stream 
direction in the potential flow. If @ is the angle between 
the x-axis and the direction of the potential stream at 


any point in the «y-surface—.e., 


TABLE 1 

Profile Slopes at the Surface 
8 m g'(O) h'(0) g*"(0 aw 
—().1988 —(). 0904 2.70 0.48 2.70 0.48 
—(. 14 —(). 0654 2.22 0.66 1.98 U. 46 
0) 0 2 02 0.76 1.55 0.29 
0.2 0.11 1.91 0.82 1.22 0. 1s 
0.5 0.33 1.83 0.84 0.90 —0.09 
- 1) 
1.0 1.0 1.76 0.81 0.53 —U.a0 
1.6 10 1.74 0.70 0.22 —0.d2 
90 . 1 —| OY 


69 0.60 O 
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Fic. 2. Boundary-layer profiles of the spanwise flow due to 
rotation 


6 = tg! u1/m (3.14) 
then the new velocities at that point are 


u* = ucosé+vsin 6 . 


* 


v* = ycosé@ — usin Of 


In the domain under consideration, however, 7 and v 
are much smaller than “ and u, respectively, so that, 
to an approximation which is consistent with Fogarty’s 
basic simplification, 


ew v™ =v — (v)/m)U (3.16) 


The new point of view does not affect the chordwise 
profile which, to the approximation involved, is identi- 
cal to the ‘“‘streamwise’”’ profile. By definition, v* is 
zero outside the boundary layer. 


From Eqs. (3.1), (3.7), and (3.8), v* is written as 





v* = [we/(m + 1)] (v/c)™*! (g — f’) - 
2w cos yx (h — f/f’) (3.17) 
Putting 
g — f’ = g*, h-—f'=h* (3.18) 
and using nondimensional quantities, Eq. (3.17) be- 
comes 
v*/wx = [&"/(m + 1)]g* — 2 cos yh* (3.19) 
or 
v*/wo = [&"+!/(m + 1)]g* — 2 cos yih* (3.20) 
The functions g* and h* are presented in Fig. 3. Table 


| also contains g*’(0) and h*’(0). The component of 
the shearing stress perpendicular to the resultant po- 
tential-flow direction can be computed from these 


quantities. 


(4) DISCUSSION OF RESULTS 


The g-profiles, Fig. 2, represent the velocity distri- 
bution for a part of the spanwise flow which is always 
directed outward. 
the resultant flow direction, also represents an outflow 


The same part, when viewed from 
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for all values of 8, as shown by the g*-profiles, Fig. 3. 
The maximum value of g* is largest for the ‘‘separation”’ 
value of 6 = —0O.1988, and gradually diminishes as 
the flow becomes more accelerated. In the limiting 
case 8B = 2, g* becomes zero, as g = f/f” in this case. 
[See Eqs. (3.4) and (3.9).] 

On the other hand, with the sign convention ex- 
pressed in Eq. (3.8), positive / corresponds to an inflow 
for positive values of x cos y. Viewed from the re- 
sultant flow direction, the plot of h* in Fig. 3 shows 
outflow profiles for high acceleration and inflow for 
retarded flow. [See Eq. (3.20).] The crossover is 
near to 6B = 0.5. 

The resultant profiles show a complicated behavior; 
not only are the g*- and /*-profiles superimposed with 
variable weight along the chord, but also the factor 
cos y has to be taken into account. As an example, 
the case cos y = | for all @ has been considered. In 
Fig. 4, the value of (v*/wC)n; is plotted against & 


for 0 < & < 1, as given by Eq. (3.20). 


Fs The point 


& = 1 has the significance that, at that point, the 
chordwise velocity “ reaches (by different power laws) 
the same value u, = wy for all m. [See Eq. (1.4).] 
With this normalization, the greatest outflow is found 
for retarded flow; the smallest values occur for 6 
around 0.2 to 0.5, while extremely strong acceleration 


* 


increases the outflow again. The value of v* is non- 
dimensionalized with wc; a comparison with 1; is of little 
value, as the spanwise boundary layer is independent 
of y, while “4 grows proportional to y. 


It may be recalled that the present example is carried 


out for the special case C = O0—1.e., v = O at the nose 
& = 0. By a shift of the axis of rotation a nonvanish- 
ing value of v appears for £ = O, and the lines of Fig. 4 


are shifted correspondingly. 

An application of the results of this paper is given 
by the use of the profiles calculated here to the con- 
struction of a momentum method for arbitrary pressure 
distributions. Some results of this investigation are 
discussed in reference 11. 
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Boundary-layer profiles perpendicular to the free- 
stream direction 
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away in the chordwise direction. If is sutticienth 







t2 ; . : large, this is possible without a noticeable increase j; 

the chordwise boundary-layer thickness. Ultimate) ( 
however, for large chords, 6 will grow stronger which 
z 08 implies increased values of I’, according to Eq. (5.3 
,£ aad | | | Thus, while the present treatment appears to be ac 
3 quate for configurations as found on propellers an¢ 
PS helicopter rotors, it can be very misleading for shor 


blades of considerable chord. 
The inclusion of the spanwise pressure gradient int 


(v 
O 
> 


the present consideration will not change the result ex 
0 oa | “0: pressed in Eq. (5.5), as the spanwise flow due to th 
O O02 04 06 O8 10 pressure gradient is of the same order, according to Ey 
& ct FE}, An 


Fic. 4. The maximum total outflow velocity in the boundary I inally, the previous rough estimate of the outflow 





~ layer perpendicular to the free-stream direction along the will be applied to some considerations of the turbulent shock 
shord. . . : fth 
ain boundary layer. In view of our incomplete knowledge . 
eee of the turbulent shearing stress, let it be assumed that 

(5) CONCLUSIONS . ; ; ; 
the spanwise and chordwise components of the shear » the 
It has been noted already by Fogarty’ that the span- are connected by a formula 

wise component of the laminar boundary layer is found Vy . 

. . . bs 2 rT Ty foes Tr\ } 0.0 . \ 
to be surprisingly small on rotating blades. The pres- 
ent analysis fully confirms this conclusion; the fol- The exponent will be discussed presently. If, again * 

- - <j ropdic 
lowing comments are intended to clarify the physical the boundary-layer thickness is determined by the f heli 
significance of this result. chordwise flow development, then ised 

The order of magnitude discussion of Fogarty” led to F 1 a 

° - 5 . 0 ™ CT;/ pl” dA 
the significant result (also announced by Prandtl!*) 
that if the chord under consideration is small compared Eqs. (5.6) and (5.7), introduced in Eq. (5.1) give 
to the distance y from the axis of rotation, then the ree i - 
‘ , ° mm tC / 9° ) 

chordwise boundary layer can be calculated from a : 
two-dimensional “‘strip’’ theory. It will be shown that The value of the exponent » is naturally | in the laminar 
this conclusion also implies the order of magnitude of the case; one might speculate that the same value is also 
velocities in the spanwise boundary layer. appropriate for a fully developed turbulent boundary 

Consider the volume of flow in the boundary layer layer. If the resultant shear is established essentiall 
with thickness 6 over a small surface portion S of the by the chordwise flow velocity, then » = 1 implies that 
blade. The centrifugal force acting on this part of the chordwise shear component is proportional to the Dimen 
the layer will be of the order pSé(U*/y). For a rough chordwise velocity component. Note that 7 > | means 
estimate of the outflow, consider this force balanced by increased outflow. The uncertainty in the value of 1 
the shearing stress in y-direction: corresponds to the same difficulty encountered for 

. “. 779 = turbulent boundary layers over yawed cylinders.*:' \ 
Sty ~ pSé6 U?/y (5.1) age 5 Sera : ; 
It is interesting to note that the boundary-layer x 

Now put Ty ~ wl (5.2) thickness for a fully developed turbulent layer varies Sub 

Pox . : with the span approximately as y~'®), in contrast 

where I’ is the velocity of the spanwise flow to be — PI EN, i a 
. to the variation with y~'? for the laminar case. This 
determined: po ae cals ee ; 
means a strong “redistribution’’ after transition. It , 
V ~ U? 8/vy (9.3) is safe to say that the strongest effects of rotation will | 
4 r ° . st 
: . : be felt in the transition region, between the laminar 
If the boundary-layer thickness is determined by the ; : ae é ; ah dimer 
Moret vier . separation point and the beginning of the fully estab- 
strip theory’ of the chordwise flow, then, at a chord é | these 
‘ lished turbulent boundary layer. 
’ x ss ; Case | 
/ | 
5 ~~ VV wc/U 5.4) betwe 
’ vel (9.4 REFERENCES 
which gives with Eq. (5.3) l Sears, W. R., Potential Flow Around a Rotating Cylindri Rec 
ryyr aha Blade, Readers’ Forum, Journal of the Aeronautical Sciences * TI 
V/U ~c/y (5.5) : Nee agaies ns . ot 
. Vol. 17, No. 3, p. 183, March, 1950 search 
pl ce ; : . : / ee > The — rey on a Rotating | + 
This is the order of magnitude of the flow velocity in Fogarty, L. E., The Laminar Boundary Layer on a an 7 " 
. 7 he Blade, Journal of the Aeronautical Sciences, Vol. 18, No. 4, PP) — Engine 
the spanwise boundary layer which has been found eee ae Se Ont 
, : 247-252, April, 1951. : nt 
before. 3 Falkner, V. M., and Skan, S. W., Some Approximate Sol The 
. . ° ° - ~ wae : ou & es 
The approximation involved in Eq. (5.4) may be tions of the Boundary Layer Equations, British R&M No. 1314, | G.N 
visualized as follows: with 6 ~ y~'”, the outflow 1930 the 
Iring 


diminishes, which is only possible if the excess is carried (Continued on page 1006) 
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An Experimental Study of the 
One-Dimensional Refraction of a Raretaction 
Wave at a Contact Surface’ 


I. J. BILLINGTON?t 


Unwersity of Toronto 


SUMMARY 


\n experimental investigation is described of the one-dimen- 
jonal interaction of a rarefaction wave and a contact surface in a 
ck tube. The rarefaction wave generated by the bursting 
he shock-tube diaphragm is utilized as the incident wave in the 





nteraction with an air and helium or air and argon contact sur- 
e. The investigation shows that the initial rarefaction wave 
the shock tube deviates considerably from that predicted by 
eal shock-tube theory. When this deviation is taken into ac- 
, however, the results of the rarefaction-wave refraction 
w good agreement with the theory. It is verified that, as 


s been 


predicted in a previous theoretical treatment of the 
problem, a compression wave is reflected from the contact surface 
if helium is used and a rarefaction wave is reflected if argon is 


ised 


SYMBOLS 


= velocity of sound 
internal energy (C,/ 
chamber length 

= pressure 

= time 

= particle velocity 

= distance 

= ratio of specific heats 

= density 


Dimensionless Ratios 


f = VW = it a 
= pi/t T;; = Ti/17 
= T/T a = (t+ 1M -1 
3 = (vy — 19/29 ii = pi/p 
Ai; = pitt; / pj T = t/L 
X = x/l 


Subscripts refer to flow regions as numbered in Fig. 1. 


(1) INTRODUCTION 


y THE PAST, the shock tube has been used for the 

study of a number of interactions between one- 
dimensional wave elements. A discussion of several of 
these interactions may be found in reference 1. The 
case where a moving wave impinges upon an interface 
between two gases of different internal energy has been 
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called refraction because, in the distance-time (x, /) 
plane, the transmitted wave shows a change in direction. 


As in the case of optical refraction, there are incident, 


transmitted, and reflected waves. An experimental 
study of shock-wave refraction has already been re- 
ported. 

The case of rarefaction-wave refraction has been 


treated theoretically in considerable detail in reference 
2 where it was shown that two results were possible, 
depending upon the ratio of the internal energies of the 
two gases separated by the interface. If this energy 
ratio -;; was greater than a critical value E..it., a com- 
pression wave was reflected and the transmitted rare- 
If, on the 
an expansion wave 


faction was weaker than the incident one. 
other hand, F;; was less than F,,;, 
was reflected and the transmitted wave was stronger 
than the incident one. The value of £,,;,, was shown 
to be given by 


| (D5 Bs) [Q —_ P3P (1 = P34" (1) 
where /;; is the pressure ratio across the incident rare- 
faction wave. The theory of reference 2 further pre- 
dicts the following relationships between P3;, /5,, and 


the reflected wave pressure ratio P;;: 





Dt 
4 


wt 
Ot 


(2) | r ee 
(1) A (4) x 


t) plane for refraction 





Fic. 1. Wave configuration in the (x, 
of a rarefaction wave at a contact surface. The bursting dia- 
phragm at the origin generates the shock wave (S), contact 
surface (C), and rarefaction wave (R The rarefaction wave 
impinges upon the interface (I) generating a transmitted wave 
(T) and a reflected wave (A). The steady states between the 
waves are numbered 
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[ V Bs(1 — P33) 3 
y] a { — (ro — 1) 


Ess : Vv ~~. + ] (2 
= (P73Ps4)” = 9 

: Bs { (P34P73)"* — 2P3s°* + 1 * 

»— 


Itn4 = Bs = (P7r3P 33)" 


Eq. (2) applies to the case of a reflected shock wave and 
Eq. (3) to the case of a reflected expansion wave. 

In the experimental investigation reported here, the 
initial rarefaction wave generated by the bursting 
diaphragm in the shock tube was used as the incident 
wave in the interaction. This wave was allowed to 1m- 
pinge upon an interface between two different gases 
initially at the same temperature. The gas combina- 
tions used were air-argon and air-helium with the initial 
wave being in air in each case. According to Eq. (1), 
the former combination should produce a reflected ex- 
pansion and the latter a reflected compression at room 
temperatures. Pressure measurements were made 
across the various wave elements to check Eqs. (2) and 
3), and measurements of stagnation temperature, mass 
flow, and particle velocity were made wherever possible 
using drum-camera and hot-wire anemometry tech- 


niques. 
(2) THEORETICAL CONSIDERATIONS 


The theoretical wave system in the shock tube is 
shown in the (x, ¢) plane in Fig. 1. The shock-tube 
diaphragm is located at x = 0 and is ruptured at time 

0 sending a shock wave §S into the channel of the 
tube followed by a contact surface C. At the same 
time, a rarefaction wave R is initiated and travels into 
the chamber of the tube. The rarefaction wave im- 
pinges upon the interface I and sets it in motion as 
shown. The transmitted rarefaction wave T continues 
on into the other gas, and the reflected wave A travels 
back through the air. As mentioned above, the re- 
flected wave may be either a compression or an expan- 
sion according to the energy ratio across the initial 
interface. The steady flow regions between the various 
waves are numbered as shown in Fig. 1. 

In idealized shock-tube theory, it is assumed that 
the initial rarefaction wave is centered at the origin so 
that all the characteristic or Mach lines meet at a point 
as shown in Fig. 1 and that the expansion is isentropic 
and one-dimensional. It is further assumed that the 
gas in region 3, expanded through the rarefaction wave, 
has the same pressure and particle velocity as the gas 
in region 2, which has been compressed through the 
shock wave. On the basis of these assumptions, it is 
possible to determine the pressure ratio P3; across the 
rarefaction wave and the values of the various flow 
parameters in region 3 as a function of the diaphragm 
pressure ratio Py. Prior to the present work, only 
limited studies had been made of the nature of the 
rarefaction wave generated by the bursting diaphragm 
so that it was difficult to assess the validity of the above 
assumptions. Thus it was felt necessary to devote 
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some time during the present investigation to the study 
of the primary wave. 
It is shown‘ that the slope of a characteristic in the 


(x, f) plane is given by 
dx/dt “+a I 


If transformation is made to a dimensionless (x, 7 
plane, where x x/L, r a;t/L, and L is the length 
of the chamber, the equation of any characteristic can 
be obtained from Eq. (4) in the following form where 


the constant is different for each characteristic. 
x = (1 a4) + (@/d4) |r + constant ” 


For an isentropic centered rarefaction wave, the con 
stant of Eq. (5) is zero for every characteristic, and the 
values of w/a, and a/ays on the last characteristic are 
uniquely determined by its slope (x/a,f).4 

When the flow conditions in region 3 behind the 
incident wave and the pressure ratio across the wave 
are known, expected conditions in regions 6 and 7 can 
be calculated. For example, P?;; can be calculated 
from Eqs. (2) or (3), and then the density and temper- 
ature ratios across the reflected wave can be obtained 
by applying the isentropic relations. The particle 
velocity in region 7 can be obtained by adding the par- 
ticle velocity in region 3 to the additional velocity in- 
duced by the reflected wave. 

In each flow region it was possible to make direct ex- 
perimental measurements of pressure, mass flow (pz), 
and stagnation temperature. These quantities are 
sufficient to determine all other flow properties in the 
region. For example, in the case of the initial rarefac 
tion wave, the equation of state can be written in the 
following form: 

P34 7 M347 34 (6) 
In region 3 behind the initial wave, the adiabatic form 
of the steady flow energy equation can be written as 


Tos T; =1+ [(v4 — 1) 2|M;? (7) 


Eqs. (6) and (7) may be solved to obtain an expression 
for the density ratio I's, in terms of the measured quan- 


tities: 
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T'34 = [P34 = V Po os 2T o3a(¥ — 1) (Az34) 7] 2T oss (S 


Similarly, the temperature, particle velocity, and other 
flow parameters can be calculated from the experi- 
mental measurements. The isentropic nature of the 
wave can be verified by applying the equation 


Py = 34"! (9) 


to the density ratio determined from Eq. (S) and com- 
paring the resulting pressure with the measured pres 


sure. 


(3) EXPERIMENTAL PROCEDURE 


The experimental work was carried out in the 5- by 
3-in. wave interaction shock tube at the Institute of 
Aerophysics. This tube is described in reference 4. 
The chamber of the tube could be divided into two 
compartments by a sliding metal shutter which could 
be quickly removed immediately before the diaphragm 
was ruptured. In the preparation for an experimental 
run, air at the desired pressure was put into the section 
between the diaphragm and the slide, and either argon 
or helium at the same pressure was put into the section 
on the other side of the slide. Withdrawal of the slide 
formed an interface between the two gases. It was 
felt that the possible eddies formed by the withdrawal 
of the slide and diffusion of the two gases before the 
arrival of the rarefaction wave were effects which did 
not appreciably change the overall nature of the inter- 
action. 

Instrumentation included a schlieren drum-camera 
system which could be used to obtain (x, f)-plane 
schlieren photographs of the interaction, a piezoelectric 
pressure gage for measuring pressure jumps across the 
various waves, and a hot-wire anemometer with which 
stagnation-temperature and mass-flow readings could 
be obtained in the various flow region. The output 
signals of the hot wire and piezoelectric pressure gage 
were recorded photographically from the screen of a 
cathode ray oscilloscope. Since the hot wire is sensitive 
to both mass-flow and stagnation-temperature changes 
which occur simultaneously in the shock tube, it was 
necessary to make a pair of shock-tube runs with the 
same initial conditions and the same wire but with 
different wire heating currents. The two output signals 
could then be used to form a pair of equations which 
for the two unknowns—stagnation 
Since wire breakage was a 


could be solved 
temperature and mass flow. 
serious problem, an operating procedure was adopted 
which was designed to extract the maximum informa- 
tion from each wire. The first run was made with the 
wire heating current very small so that the wire was 
sensitive only to temperature changes.’ If the wire 
broke after the first run, the assumption that the wire 
current was approximately equal to the stagnation 
temperature of the flow could be used to make an ap- 
proximate determination of the flow temperature. If 
two successful runs were completed, the stagnation 
temperature and mass flow could be determined by the 
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exact solution of the two equations mentioned aboy, 
Hot-wire and piezo-gage measurements were taken j; 
regions 3, 6, and 7 (see Fig. 1) for both gas combinations 
and at various diaphragm pressure ratios in the rang 
i< Pa < @. 

In addition to the investigation of the various floy 
regions, the hot wire and pressure gage were used t 
study the incident wave itself. An extended time bas 
was used, and profiles were calculated from the oscillo 
scope traces. The shape and duration of these profiles 
could be compared with those for an ideal centered 
rarefaction wave. 

The reflected wave, in which the density changes ar 
small, was not apparent on the drum-camera records 
The front of the contact region separating regions 6 and 
7 did appear, and its velocity could be measured from 
the films. Since the velocity of this front is assumed to 
be identical with the particle velocities in regions 6 and 
7, the drum-camera measurements could be compared 
with the values calculated on purely theoretical grounds 

As the diaphragm-pressure ratio is increased, the ex 
pansion waves in the chamber become more and mor 
extended and the duration of the steady states between 
the waves becomes correspondingly shorter. This 
factor was found to place an upper limit on the strength 
of the waves which could be studied. The channel and 
chamber lengths and the positions of the slide, hot wire 
and pressure crystal had to be adjusted according to 
the wave strengths and the measurements being made 
Typical tube configurations are shown in Fig. 2. 

Fig. 3 shows a block diagram of the electronic equip- 
ment with the hot wire and piezoelectric pressure gage 
each connected to a beam of a twin-beam oscilloscope 
The hot-wire anemometer is of the constant current 
type with a compensating amplifier. Signal gener- 
ators are provided for the calibration of the amplifier 
and of individual hot wires. The oscilloscope was 
allowed to sweep only once for a given run and was con- 
trolled by a second piezoelectric crystal placed down- 
stream of the diaphragm and activated by the initial 
shock wave. The sweep triggering pulse could be 
passed through a time delay unit as shown. 

A detailed discussion of the apparatus and the ex 
perimental techniques employed may be found in refer- 
ence 3. A composite microscope photograph of 
Wollaston-process hot wire magnified approximately 
50 times is shown in Fig. +. The 0.0001-in. diameter 
platinum wire is obtained with a large silver sleeve and 
is soldered in a semicircular form between the tips ol 
two sewing needles as shown. The silver is etched off a 
portion of the loop leaving the platinum wire bare. 

The use of the hot-wire anemometer in shock-tube 
studies was first suggested by Kovasznay and Do 
sanjh'! '? who developed experimental techniques and 


methods for the reduction of data. 


(4) NATURE OF THE INITIAL RAREFACTION WAVE 


It became apparent early in the investigation that 
the pressure, stagnation temperature, and mass flow 10 
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the region behind the initial rarefaction wave deviated 
considerably from the theoretically predicted values 
except for very low diaphragm-pressure ratios. Figs. 5, 
6, and 7 show the measured values for the pressure, 
stagnation-temperature, and mass-flow ratios, respec- 
tively, as a function of diaphragm-pressure ratio. Each 
experimental point on these graphs represents an 
average of a number of experimental runs at the same 
pressure ratio. For comparison, each graph also shows 
a theoretical curve computed on the basis of an ideal 
isentropic centered rarefaction wave. 

The pressure ratio across the incident wave was 
found to be larger (see Fig. 5) than that predicted by 
This indicated that the wave was weaker than 
expected. The deviation is negligible for low dia- 
phragm-pressure ratios but becomes very marked for 
The stagnation-temperature 


theory. 


the higher values of P4). 
ratios, on the other hand, are somewhat lower than pre- 
dicted by theory as can be seen in Fig. 6, although the 
discrepancy is not large anywhere in the range investi- 
gated. The accuracy of the mass-flow measurements 
was felt to be less reliable than the pressure and tem- 
perature measurements because of a functional depend- 
ence of the mass flow upon a difference between two 
small measured quantities of the same order of magni- 
tude. This was reflected in a relatively large scatter 
in the measured mass-flow values. The average values 
of the mass-flow ratio (Fig. 7) indicate a considerable 
discrepancy between theoretical and measured condi- 
tions, especially at higher diaphragm-pressure ratios. 

The deviations cannot be explained on the assumption 
that the actual waves are one-dimensional isentropic 
centered rarefaction waves but merely weaker than 
expected for a given diaphragm-pressure ratio. If this 
were the case, the discrepancies on Figs. 5, 6, and 7 
would be consistent. In fact, they are not consistent 
as can be seen from consideration of the case where the 
According to the 


diaphragm-pressure ratio Py, = 3.5. 
pressure measurements (Fig. 5), this wave is somewhat 
weaker than expected. According to the stagnation- 
temperature measurements, however, (Fig. 6), this 
wave is considerably stronger than expected. 

Using the measured values shown in Figs. 5, 6, and 7, 
the density ratio across the incident wave was deter- 
mined from Eq. (8) for each diaphragm-pressure ratio. * 
These calculated density ratios were found to be in good 
agreement with direct density measurements made with 
a chrono-interferometer by Mack. When Eq. (9) was 
applied to these density ratios, it was found that the 
resulting pressure ratios were in very good agreement 
with the measured values of Fig. 5, indicating that the 
initial rarefaction wave is in fact isentropic. From 
the above measurements, the remaining flow param- 
eters such as particle velocity and static temperature 
were calculated in region 3. The results of these cal- 
culations are illustrated graphically and discussed in 

* These calculated density ratios are shown in Fig. 23. A 
theoretical curve computed on the basis of an ideal isentropic 
centered rarefaction wave is also shown for comparison. 
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some detail in reference 3. The particle velocities ca]. 
culated from the measurements in region 3 were found 
to deviate from those calculated from the one-dimep. 


sional characteristic equation 
~ = [2a ‘(¥4 — 1) ] = = (2a, je = l (10 


It was felt that this discrepancy was not due to any 
systematic error in the measurements or to the use of 
the wrong heat-loss equation for the hot wires becaus: 
the particle velocity could be obtained from the meas 
ured quantities by two independent methods. Either 
the measured mass flow from the hot-wire records 
could be divided by the density ratio calculated from 
Eq. (8), or the measured stagnation temperature could 
be divided by the static temperature, calculated from 
the measured pressure according to isentropic theory, 
to yield the Mach Number according to Eq. (7) and 
hence the particle velocity. The latter method is inde 
pendent of the measured mass-flow values which were 
considered to be the least accurate of the measurements 
taken. 

A dimensionless (x, 7)-plane diagram of the initial 
rarefaction wave was constructed for the case Ps, 
2.0 and is shown in Fig. 8. On the assumption that 
the first characteristic passes through the origin, the 
value of 7 on the first characteristic at any station can 
be determined from Eq. (5) since u 0 anda = a. 
Then, from the measured duration of the wave at that 
station as determined from the extended time base 
studies with the hot wire and pressure gage, the value 
of 7 for the last characteristic at the same station can be 
determined. Using the known values of w/a, and 
a/a; at the tail of the wave, the slope of the last char- 
acteristic can be calculated according to Eq. (5). Thus 
a value of 7 and a slope are known at a given station (x), 
and the last characteristic can be drawn as shown in 
Fig. 8. The last characteristic of the theoretical cen- 
tered wave for the same diaphragm-pressure ratio 1s 
also shown for comparison. It can be seen that the 
two waves have considerably different forms and that 
the real wave is not centered at the origin. 


(5) RAREFACTION-WAVE REFRACTION 


A quantitative examination of the photographic 
records shows that the general wave system set up in 
the shock tube as a result of the rarefaction-wave- 
contact-surface interaction is of the type predicted by 
the theory. In the air-helium case, a compression wave 
is reflected as shown in Fig. 9, and, in the air-argon case, 
an expansion wave is reflected as shown in Fig. 10. 

Fig. 9 is a typical piezoelectric pressure gage record 
in the air-helium case showing the incident rarefaction 
wave (R) followed by the reflected compression wave 
(S’). The moving contact surface does not show up 
on the piezo-gage trace. The crystal station is indi- 
cated on the accompanying sketch of the (x, ¢) plane. 
In the air-argon hot-wire trace in Fig. 10, the incident 
rarefaction wave (R) is seen to be followed by the re 
flected rarefaction (R’) and the moving contact 
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Fic. 18. Variation of particle velocity after the interaction 
(Uz, = Uses) with the pressure ratio across the incident wave 


(P34). Case: air-helium. The experimental points are from 
drum-camera records, and the calculated points are from hot- 
wire and piezo-gage measurements 
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Fic. 19. Variation of the stagnation temperature behind the 


reflected compression wave (7073) with the stagnation-tempera- 
ture ratio across the incident wave (7 3;). Case: air-helium 
(The experimental points determined from one hot-wire run are 
shown as triangles, while those determined from a successful pair 
of runs are shown as circles. ) 
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Fic. 22. Variation of the mass-flow ratio across the trans- 
mitted wave (Ags) with the stagnation-temperature ratio (7 34) 
across the incident wave. Case: air-helium. 
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transmitted wave (7 4) with the stagnation-temperature ratj 
( The experi- 
mental points determined from one hot-wire run are shown as 
triangles, while those determined from a pair of successful runs 


across the incident wave (73s Case: air-helium 


ire shown as circles. ) 
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Fic. 21. Variation of the mass flow behind the reflected com 


pression wave (A74) with the stagnation-temperature ratio across 


the incident wave (7%;;).. Case: air-helium 
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Fic. 23. Density ratio (I'3,) across the incident rarefaction 
wave. The values shown are calculated from hot-wire and 
piezoelectric pressure gage readings and are compared with 
chrono-interferometer measurements by Mack. 
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ONE-DIMENSIONAL REFRACTION OF A RAREFACTION WAVE 1005 


region (C). The reflected rarefaction wave from the 


interaction has only about one-tenth the amplitude of 


the incident wave. In both cases it can be seen that 
the steady state following the contact region is ulti- 
mately washed out by the rarefaction wave (R”) re- 


flected from the closed end of the shock-tube channel. 


The measured results are shown in Figs. 11 through 
16 in the air-argon case and in Figs. 17 through 22 in 
the air-helium case. The pressure-gage and drum- 
camera readings are plotted against incident rarefac- 
tion-wave pressure ratio, while the hot-wire readings 
are plotted against incident rarefaction-wave stagna- 
tion-temperature ratio. It can be seen from Figs. 11 
ind 17 that the measured pressure ratios agree very 
well with the theoretical values which were calculated 
The solid curves 


irom Eqs. (3) and (2), respectively. 


in the remaining figures represent the expected real 
flow values based on the actual conditions in region 3 
The dotted 


curves, which are shown for comparison, are calculated 


ud calculated from the isentropic theory. 


from the same theory but are based on the conditions 
which would exist in region 3 in the case of ideal 
shock-tube flow. Details of the computation of these 
curves are given in reference 3. Distinction is made 
between experimental points obtained with the probe 
initially downstream of the diaphragm station (in re- 
gion 4) and those obtained with the probe initially up- 
In the 


case of stagnation temperatures, distinction is also made 


stream of the diaphragm station (in region 5). 


between results obtained from a successfully completed 
pair of shock-tube runs and approximate results ob- 
tained from hot wires which broke after only one shock- 
tube run. 

It can be seen that the measured results are generally 
in good agreement with the real flow curves, particularly 
in the case of stagnation temperature. In the case of 
mass flow, considerable scatter in the experimental re- 
sults is evident, especially in Figs. 15 and 22. Figs. 
i2 and 18 show that the measured values of the contact 
iront velocity tend to be slightly higher than the real 
flow theory predicts. This discrepancy may be due 
to the fact that mixing of the two gases causes the 
thickness of the contact region to increase with time so 
that the apparent velocity of the front of the region is 
greater than the velocity of the adjoining flow regions. 


From the piezoelectric and hot-wire measurements 
in regions 6 and 7, it is possible to calculate the density 
ratios across the transmitted and reflected waves using 
relations equivalent to Eq. (8). The appropriate isen- 
tropic relationships in the form of Eq. (9) may then 
be applied to these density ratios, and the resulting 
pressure ratios can be compared to the measured piezo- 
gage values to check the isentropic nature of the waves. 
This test was applied to a number of representative 
points for both gas combinations. The average dis- 
crepancy between calculated isentropic pressure ratios 
and experimentally measured pressure ratios was 
found to be 2.2 per cent, indicating that all the wave 


elements are isentr¢ ypic. 


Using the density ratios determined as outlined above 
and the measured mass-flow values, the velocity ratios 
across the transmitted and reflected waves could be 
calculated. This was done for several representative 
points, and the results are superimposed on Figs. 12 
and 18 for comparison with the drum-camera results 
and the theoretical curves. It is evident that these 
points do not agree too well with the other values, 
probably because of the large scatter in the mass-flow 


results (see Fig. 15, for example 
(6) DiIscussION 


The initial investigations of the incident rarefaction 
wave yielded interesting results which appear to estab- 
lish the fact that the wave generated by the bursting 
diaphragm is isentropic but is not centered anywhere 
in the (x, ¢) plane. This result is substantiated in part 
by recent drum-camera runs using Ronchi grating 
schlieren which indicated’ that the rarefaction was not 
centered at the origin. The calculated density ratios 
across the incident wave compare well with the meas- 
urements of Mack® despite the fact that these measure- 
ments were made in a small diameter (1*/s in.) circular 
tube. Previous studies of the pressure profiles through 
rarefaction waves by Huber, Fitton, and Delpino® indi- 
cated agreement with ideal flow theory, but that work 
was confined to diaphragm-pressure ratios less than 
2.5, and the present work also shows good agreement 
with theory in that range. It should also be noted that 
the incident wave does not satisfy the one-dimensional 
characteristic equation. 

Although the measurements of the various ratios 
across the incident wave check well among themselves 
and with the work of other observers, they do not agree 
with measurements made‘ in that portion of region 3 
which lies in the shock-tube channel. For example, it 
has been verified that the pressure ratio 2, (see Fig. 1) 
across the incident shock wave is in good agreement 
with ideal shock-tube theory for diaphragm-pressure 
ratios such as used in this work. A study of the contact 
region’ between states 2 and 3 showed that no pressure 
discontinuity existed there. However, the present 
measurements in the shock-tube chamber (Fig. 5 
indicate deviations from the ideal flow theory. The 
main purpose of the present measurements, however, 
was the establishment of a basis for an understanding 
of the refraction problem, and the effects at the dia- 
phragm station will therefore require further study. 

The results of the measurements in the various flow 
regions generated by the interaction of a rarefaction 
wave with a contact surface appear to be in good agree- 
ment with the theory when account is taken of the 
actual form of the incident wave. The numerical re- 
sults indicate agreement with the theory within the 
limits of the experimental error for the two gas combi- 
nations investigated. It appears therefore that the 
analysis of reference 2 can be relied upon to give ac- 
curate results in other cases also. It should be noted 
that relatively large variations in velocity in region 3 
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On Turbulent Flow Near a Wall 


E. R. van DRIEST* 


North American 


SUMMARY 


A theory which provides a continuous velocity and shear 
jstribution for turbulent flow near a smooth wall is developed 
The analysis also forms the basis for the theoretical calculation 
f the velocity profiles and resistance owing to roughness or vortex 
generation The theory checks well with experimental data 


SYMBOLS 


distance from wall 
pipe radius 
D pipe diameter 
i mean local velocity parallel to wall 
= velocity fluctuations parallel and normal to flow 
mass density 
coefficient of viscosity 
r = shear stress 
velocity correlation coefficient 
mixing length 
= universal constant in/ = Ay 
modified universal constant 
eddy viscosity 
= size of roughness 
= friction factor =8r,,/pV? 
average section velocity 
Reynolds Number = VD/1 
= kinematic viscosity 
= subscript denoting wall 
= subscript denoting turbulence 
= subscript denoting turbulence 
= subscript denoting maximum 


= UuVr p 


iol 
€¢ £:°& 
me | 


= V 7,,/p°a/v 


constant 


INTRODUCTION 


, ‘HE PURPOSE of this report is to present an analysis 
which will yield a continuous velocity and shear 
smooth wall. 


for turbulent flow 


The results of such an analysis should be particularly 


listribution near a 
useful in the study of convective heat transfer by 
liquids over a wide range of Prandtl Number. The 
report pertains to incompressible fluids only. 

It is found that the theory fits the measured velocity- 
profile data for smooth walls. The analysis also leads 
to the theoretical calculation of velocity profiles and 
iriction factors for rough pipes and ‘or vortex generators. 


ANALYSIS 


When the Navier-Stokes equations are written in 


ter ~ ” “4° ° e 
‘fms of mean velocities and fluctuations from the 
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y Lviation, Ine. 


mean and then averaged with respect to time, rearrange- 
ment of the resulting equation allows the total shear 
stress 7 to be identified as 


7 = p(On/Oy) — pu’v’ (1) 


where # is the mean velocity parallel to the wall, u’ 
the instantaneous fluctuation of velocity in the direc- 
tion of #, v’ the cross fluctuation of velocity in the 
direction normal to the wall, y the length scale normal 
to the wall and measured positive from the wall, p 
the density of the fluid, and u the coefficient of viscosity 
of the fluid. 


] 


Eq. (1 


The first term on the right-hand side of 
represents the effect of viscosity on the mean 


flow whereas the second term is the Reynolds stress. 
In turbulent flow away from a wall, the Reynolds 
stress is of considerably greater magnitude than the 
the wall is 


approached, the greater becomes the role of the viscous 


viscous stress; however, more a smooth 
stress until finally, at the wall, viscosity predominates. 

The viscous effect of the proximity of the wall may 
be estimated in the following manner: 

Consider an infinite flat plate undergoing simple 
harmonic oscillation parallel to the plate in an infinite 
fluid. As was shown by Stokes,' the amplitude of the 
motion diminishes from the wall as a consequence of 
the factor exp(—y/A) where A is a constant depending 
upon the frequency of oscillation of the plate and the 
kinematic viscosity v of the fluid. Hence, when the 
plate is fixed and the fluid oscillates relative to the plate, 
the factor [1 — exp(—y/A)] must be applied to the 
fluid oscillation to obtain the damping effect on the 
wall. Now, Prandtl, the total 


stress for turbulent flow is written as 


according to shear 


; 


rT = p(On/Oy) + prVu"?-Vv" (2) 


where 7 is the correlation coefficient defined by u’v’ = 
—rVu"-V 9". Introduction of the mixing lengths 
/, and /, defined by Vu’? = /, On/O0v and Vv” = 
/, Ou/Oy, respectively, then gives 


u(On/Ov) + pril.(Onm Oy)? (3 


T= 


The quantity r/,/, is finally lumped into one overall 
length / to yield 


rT = p(Ou/Oy) + pl*(Ou/dy)? (4) 


Furthermore, when it is assumed that / = Ay, where 


K is called the universal mixing constant, Eq. (4 


becomes 
rT = p(On/Ov) + pK*y*(0n/dy)? (5) 


which is supposed to represent mean fully developed 


turbulent flow near a wall. However, such fully 
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developed turbulent motion occurs only beyond a According to Eq. (1), the Reynolds str 
distance sufficiently remote from the wall that the very obtained from 
eddies themselves are not damped in turn by the 
nearness of the wall. Indeed, near a wall, the damping SS eee Th 
factor would be [1 — exp(—y/A)] for each mean in which 7, = —pu’v’. Hence, with Eq. (9 
velocity fluctuation in Eq. (2), and hence Eq. (5) 

T, = T — Ty(OUx/OVx ) 


should be modified to become 


whence Ti/Te = (17/Tr) 


r = w(On/Oy) + pK*y*[1 — exp(—y/A)]?(Ou/Oy)? (6) 


— (Oux /OVs 


The eddy viscosity ¢€ is obtained from 


in order to take into account the mean motion all the 


way to a smooth wall. r = uw(On/Oy) + €(OU/Oy) = 


SS T 18 


(MT €)(Tw/ Mh) * (Olt / Oy, 






One could argue that the presence of the wall modifies r 
the universal constant so that 
; rm so that e/a = [(7/ty)/(Oux/OVe) |] — 17 
x = K[l — exp(—y/A)] (7) 
or that the mixing length must be changed to Smooth Wall 
1 = Ky[l — exp(—y/A)] (8) For boundary-layer flow with zero pressure gradient 
: ; : aay : Or/Oy = 0 at the wall and therefore 7 & 7,, near the 
It is convenient to write Eq. (6) in dimensionless . : , - * 
: oe ; wall. Hence, Eq. (10) yields 
form. Thus, putting 
es / / OUx 2 
Ux = U/V Tw/ p Va = V Tw/p°y/v (9) = ; 18 
; Ove 1+ V1 + 4K°%yx2[1 — exp(—ye/Ax)] 
where 7, is the shear at the wall, Eq. (6) becomes 
oe A whereupon 
T/T. = (OUx/OVx) + K?yx?[1 — exp(—ys/Ax)]? X 
9 P 2d Vx 
(Oux/OVx)? (10) gt ms | ; : 
3 : J . . v0 I+Vl1ls+ 1K? yx7[1 — Can — Va A, 
in which Ax is a constant of the turbulence. Further- | 
14 
more, 
k = K[l — exp(—y«/Ax)] (11) Then, from Eqs. (15) and (17), respectively, 
ly = Kyx[1 — exp(—yx/Ax)] (12) T,/Te = 1 — (Oux/OV«) 20) 
where = W/7,/p:l/v. e/u = [1/(Oux/Ovx)] — 1 2] 
30 r 2 
TURBULENT 
(PARTLY ROUGH ) 
+ EO. 33 


TURBULENT ky 
(SMOOTH WALL ) 


LAMINAR EQ 19,28 


20 Ux= Vg -4 + - 
Vy] END OF WALL 
INFLUENCE 
/ Yy=60 . b 
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Fic. 1. Semilogarithmic plot of velocity profiles for turbulent flow near smooth and rough walls 
periment for a smooth pipe 
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TURBULENT 
For large vx, Eq. (19) becomes 
ux = const. + (1/K) In y, (22) 
which is the von Karman law for fully turbulent flow 
i.e., where viscosity does not affect even the mean flow. 
The above analysis applies to a smooth wall under 
which condition the proximity of the wall has a stabiliz- 
ing effect upon the eddies owing to viscosity. Of 
course, the opposite could be argued—namely, that 
the wall loses its control of the fluid elements with 


distance from the wall according to the same law. 


Rough Wall (Vortex Generation) 


[The stabilizing effect of the wall can, however, be 
nullified by introducing an artificial mixer at the wall 
When sufficient roughness is intro- 
near the wall (vortex 

A,) should disappear 


viz., roughness. 
duced to stir up the motion 
generation), the factor exp(— yx 


fom the above equations so that Eqs. (11), (12), 
and (19) should become, respectively, 

K K (23) 

K yx (24) 

Oux / OV 2/{1 + V1 + (2Ky%)?] (25) 


which integrates to 


L fl — V1 + (2K ye)? 
K 2K ys 


In | 2K, +Vi1+ Ky’) | (26) 


For large yx, Eq. (26) becomes 


ux = const. + (1/K) In yx 


which is again the von Karman law for fully developed 
turbulent flow. 

It is expected that Eq. (26) should represent the 
beginning of fully turbulent flow clear to the wall 
1.¢., flow which feels only the effect of viscosity on the 
mean motion and not a viscous effect of the nearness 
of the wall on the individual eddies. 

The above results for flow in the vicinity of both 
smooth and rough walls are also approximately valid 
with a pressure gradient in the direction of the flow 
because the shear stress near the wall is approximately 
equal to the wall stress. 


COMPARISON WITH EXPERIMENT 


Velocity Profiles 


Fig. 1 shows a semilogarithmic plot of the mean- 
velocity data obtained by Laufer? for fully developed 
turbulent flow in a smooth pipe at two Reynolds 
Numbers—namely 50,000 and 500,000—based on the 





diameter of the pipe. Also plotted in the figure is 
Eq. (19) when K = 0.4 and Ay = 26. It appears 
that the theory follows the data quite well 
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Fic. 2. Comparison of theory and experiment for turbulent flow 


near the wall of a smooth pipe 


For a smooth wall, the asymptot'c curves at the 
wall and out in the so-called fully developed turbulent 


flow are, respectively, 


Us Ve 
Ux 5.24 + 5.75 logy ys 

Fig. 2 is a Cartesian-coordinate plot of the Laufer 
data and Eq. (19) for the range yx < 70. 

It should be mentioned here that, in order to fit 
Eq. (19) to Nikuradse’s smooth-wall data,* it would be 
necessary that Ay = 27. 

The velocity profile for the beginning of complete 
roughness, Eq. (26), is plotted in Fig. 1 with A 0.4. 
The asymptote in the stream is 


»+) 


Ux = —1.325 + 5.75 logy vs (29) 


Eq. (26) is plotted also in Fig. 2. Of particular interest 
is the fact that only one constant, A, appears in Eq. 
(26). 

It is evident from the smooth-wall curve of Fig. | 
that the viscous damping effect of the wall extends 
out to about yx = 60. Therefore, it is expected that 
any roughness elements should also extend to about 
60 before they completely nullify the viscous 
Thus, if there are no viscosity 
VTn/p X 


Ve = 
influence of the wall. 
effects for roughnesses greater than ky 
k/vy = 60, where & is the average roughness size, 
then the general velocity profile beyond the roughness 
protuberances would be, from dimensional analysis, 


ux = const. + (1/KA) In (y/k) (30) 
= const. + (1/A) Inky + (1/K) In ye (31) 

so that, from Eq. (29) when A = 0.4 and kx = 60, 
Ux = 8.95 — 5.75 logis Rx + 5.75 logio Vx (32) 


The two regions of flow, one under viscous influence 
of the wall and the other with a wholly rough wall, 
are indicated in Fig. 1. In the region of viscous 
influence of the wall, the roughness height ky is less 
than 60—1i.e., for kx < 60, the nearness of the wall 
still shows some effect through viscous damping. 

An estimate can be made of the velocity profile 
when the roughness is within the v scous influence of 


i.e., kx < 60. It was seen above that the 


the wall 
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beginning of the complete roughness effect was felt 
on the flow when the roughness protuberances pro- 
‘ected out sufficiently far from the wall to disrupt the 
Thus, the factor exp(—y« = 
this 


be expressed mathematically by adding a local vortex- 


viscous effect of the wall. 


4A,) was artificially eliminated. However, can 


generation factor exp(—%/26) to the damping factor 


-exp(—y«/26). But the vortex-generation factor 
should grow with the size of the roughness, and, 
therefore, such a factor should have the form 
exp(—60ys 26kx) so that, at vx = ky 60, the dis- 
turbance factor owing to roughness just offsets the 


the 
The velocity profile 


damping factor. Obviously, as kx changes, 
roughness factor remains similar. 


including roughness then becomes 


[ 2dy 
i Yo 1+V1+4K"y,?[1 — exp(—ys/26) + 
exp( — 60yx /26k) | 


99 
(eo) 


so that the damped universal constant and damped 


mixing length become, respectively, 


K{1 exp(—yx/26) + exp(—60y%/26ks)] (384) 
Ky [1 — exp(—y«/26) + 
exp(—G60yx /26k, (35 


Some velocity profiles outside of the roughness elements 
are plotted in Fig. 1. It is seen that the smooth wall 
results when k, 0. On the other hand, if Eq. (32) 
is assumed to hold for all sizes of roughness, then 
that, not 

() but rather that k, 


(2S) and (32) 


it will be found for smooth effects, it is 


necessary that kx t because, 


from Eqs. 


5.24 8.95 — 5.75 logy ks 
whence ky 1. It seems more logical, however, that 
k, = O rather than k, t represents a smooth wall. 


Friction Factor 


It is readily shown that the mean and maximum 


velocities of the semilogarithmic profile are related by 


(Umar — V)/WTx/p (3/2K) 3.70 (36) 

where |’ is the average section velocity and AK 0.4. 
Also, the friction factor f for pipes is given by 

, = ST y ‘pV? (37) 


Hence, from Eqs. (28)—with y = a, (36), and (37), 


there results for a smooth wall 


V8 /f 1.49 + 5.75 logy [~/7./p- (a/v) | (38) 
1.49 + 5.75 logy [~/f/a-(V/2)-(D/v)] (39) 

1.49 + 5.75 logy (1/42) + 5.75 log; (Ry f) 

(40) 

—2.84 + 5.75 logy (Rv f) (41) 


In which D 2a and R VD/v. 
sion by 1/8, the friction law for a smooth wall becomes 


Hence, upon divi- 
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l/Vf = —1.0 + 2.04 logy (Rv f) (42 
However, the test data of Nikuradse* indicate that 

the constants should be adjusted so that 
l/-Vf —O0.8 + 2 logy (RY f (43 


For the beginning of complete roughness, Eq. (29 
gives 


9 


l/Vf —3.3 + 2 logw (Rf 
Eqs. (43 with 
curves faired through the rough-wall data of Nikuradse.* 


and (44) are plotted in Fig. 3 along 


It is seen that the semitheoretical formula, Eq. (44), 


crosses all roughness curves at approximately the 


end of transition for each relative roughness. This 


result is particularly interesting because only one 
constant, A, entered in the derivation of Eq. (44). 
Thus, the concept of generation of vortices by roughness 
appears justified. 

The general formula for complete roughness (inde- 
pendency with Reynolds Number) is now deducible 


from Eqs. (32), (36), and (37), whence 


V8/f 5.15 — 5.75 logy ke + 5.75 logy ay (45) 
where d, V T./p a/v. Division by V8 gives 
1/~/f 1.82 — 2.04 logy (ax /Rx) (46 


which agrees well with the experimental formula —viz., 


l/vV ff 1.74 — 2 logy (ax/ks (47 
representing the horizontal relative roughness lines in 
Fig. 3. It is noted that the only constant necessary 


to be evaluated experimentally in the derivation of 


Eq. (46) was A in / Ky. 
It is also possible to estimate theoretically the 
friction curves in the transition region of Fig. 3—1.e., 


curves of various relative roughness connecting the 


curves of Eqs. (48) and (44). Thus, one such curve, 
A, is obtained for a/K 60 through use of 


with arbitrary ky and the 


curve 
Eqs. (33), (36), and (37) 
relation 


( 1S) 


RV} 


Curve A has the proper trend. 


Other Properties Near a Wall 


Figs. 4 and 5 show the distributions of the damped 
universal constant and of the damped mixing length 
according to Eqs. (34) 
K 0.4. 

Fig. 6 gives the corresponding eddy viscosity ac- 
Also plotted in Fig. 6 
are the assumptions of Taylor? and Prandtl,' 
Karman.’ The shear distribution, Eq. (20), is 
shown in Fig. 7. 


and (35), respectively, for 


cording to Eq. (21) for r = Ty. 
and of 
von 
the cross 


Direct measurements of 


correlation u’v’ are apparently not available for yx 
less than about 100. 


(Continued on page 1036 








Laminar Jet Mixing of Two Compressible 
Fluids with Heat Release’ 


S. I. PAIt 


Unwersity of Maryland 


SUMMARY 


The laminar jet mixing problems with heat release have been 
formulated. A general discussion of the solution of these prob- 
lems is also given. The important parameters of these problems 
are brought out. Some specific cases of the jet mixing problem, 
such as jet mixing of one compressible fluid, isothermal jet mix 
ing of two compressible fluids, and isovel jet mixing of two com- 


pressible fluids with heat release, are discussed in detail 


(I) INTRODUCTION 


y I VHERE ARE MANY practical flow problems involving 
the jet mixing process in which mass, heat, and 
momentum transfers occur in free jet flow with or 


without heat release. The jet mixing process without 


heat release has been extensively investigated. A 
considerable amount of valuable information about 
this basic phenomenon has been obtained.' The 


jet mixing problem with heat release, however, is not so 
well understood as that without heat release. 

Most of the knowledge of jet mixing with heat re- 
lease is still on a purely empirical basis.” One of the 
most important and interesting problems which in- 
volves jet mixing with heat release is the combustion 
process in a free jet which is important in the design of 
the thermal jet propulsion system and rocket motor. 
Recently, Marble that this 
problem of combustion in jets may be investigated 
analytically by the help of boundary-layer approxima- 
tion which is essentially the usual method used in fluid 


and Adamson*® showed 


dynamics for jet mixing problems.‘ They investigated 
the ignition and combustion in a laminar mixing zone 
of two uniform streams. In this paper, we are going to 
formulate the jet mixing problems with heat release in a 
more general manner which includes the Marble-Adam- 
son problem as a special case. 

The flow of jet may be divided into two regions——the 
mixing region and the core region. Fig. 1 shows a 
typical interferogram of a two-dimensional supersonic 
jet under full expansion without heat release. In this 
figure we see clearly that there are two distinguishable 
The core region is that portion in the center 
In the core 


regions. 
of the jet and near the exit of the nozzle. 
regions viscous effects are negligible, and the fluid may 
be considered to be inviscid. The first part of the mix- 
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ing region is on the boundary of the jet near the exit 
the nozzle. This mixing region widens as the flow goe 
Far downstream, the whole jet will be 
the the 
phenomena are important. There are two types ; 
the molecular transfer 


downstreain. 


mixing region. In mixing region, transfer 


transfer phenomena ind th 
turbulent transfer. In general, these two types exist 
in tandem. In the first part of the mixing zone, th 
flow is laminar and the molecular transfer is predomi 
nant, while downstream, the flow is turbulent and turby 
lent transfer is predominant. In this paper, we con 
sider only the laminar jet mixing region of jets. 

In the flow of the mixing region of a jet, Prandtl’s 
Th 


Thus 


boundary-layer approximations are applicable. 
fundamental equations are greatly simplified. 

the mixing region of a jet is sometimes called the free 
We are going to consider the jet 
mixing of two fluids—i.e., the fluid in the jet is differ 
ent from that in which it is injected and with which it 
mixes. The essential feature of the jet mixing of tw 
fluids is that the mixture is a heterogeneous substance 
velocit} 


boundary layer. 


Besides the ordinary dependent variables 
components, pressure, density, and temperature oi 
a new variable, the concen 
Therefore, 1! 


ordinary fluid dynamics 
tration of one fluid in the other, enters. 
the fundamental equations, one has to add another re 
lation to determine this new variable—1i.e., the equa 
tion diffusion. The fundamental 
laminar jet mixing of two compressible fluids without 


of equations for 


heat release or chemical reaction have been studied 
by the author.4 Here we shall extend them to includ 
the effect of heat release due to chemical reactions 
The basic assumptions for this analysis are that th 
mixture may be considered as a continuous mediuz 
that the fluids are assumed to be perfect gases, and that 
the chemical reaction proceeds according to the Ar 
rhenius rate law and no back reaction is considered. 


(II) FUNDAMENTAL EQUATIONS 


We consider a steady two-dimensional laminar flow 
of a jet of a combustible gas mixture (gas 1) which 1s 
injected into and mixed with a uniform stream oi 
another gas (gas 2) which is the combustion product 0! 
the combustible mixture 1). The fundamental 
equations which govern the flow in the mixing regio! 


(gas 


(Fig. 2) of the gases 1 and 2 are as follows: 
(a) Equation of State 


Temperature T—The temperature T of the mixture 's 
considered as one of the dependent variables in the 
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resent problem. We assume that the temperatures of 
ises at a given point are the same and equal 
This 


consistent with our assumption that the mixture may 


the two ¢ 
to the temperature of the mixture at that point. 


he considered as a continuous medium. 


Pressure p—The pressure p of the mixture is con 
sidered as another dependent variable in the present 
problem. Since each constituent of the mixture is 


msidered as a perfect gas, we express here Dalton’s 
law as follows 

The volume of a mixture of two or more perfect 
gases is simply the sum of the volumes which they 
would occupy if each exerted the same total pressure. 
Den Let N the total 


number of molecules of the mixture of gases 1 and 2 


and Concentration be 
per unit volume under the pressure p and temperature 

N, is the corresponding number of molecules of 
is 1 per unit volume and .V, that of gas 2. We have 


N Vit iN 


We define the concentration c of gas 1 in the mix 


ture aS 


iid then the concentration of gas 2 in the mixture 1s 
1 —<¢ N./N (3) 


Let m, be the mass of a molecule of gas | and m: 
that of gas 2. The mass density p of the mixture under 


the pressure P and temperature 7 is 


p mN, + moN Nm.(Bc + 1 (+) 


where 


3 m,/Ms) —1 


We also define the mass concentration of gas | in the 


mixture as 


k m,.N;/p m,c/|m(Bc + 1 (5) 


When the molecular weights of these two gases are 
equal, then 8 = 0, c = k, and the density of the mixture 
1S independent of the concentration. In general, the 
density p of the mixture is a function of concentration. 


The equation of state of the mixture is then 
pb = NRT = pR'T (6) 
where 


R' = R/|[m2(Be + 1) gas constant for the 


mixture (7) 


The total number .V of molecules per unit volume of 
the mixture depends on the pressure p and tempera- 
ture 7’ but is independent of the concentration. 


(b) Equation of Continuity 


Since the mixture is considered as a definite fluid, the 
condition of conservation of mass gives the equation of 
continuity. For two-dimensional steady flow, we have 


the equation of continuity as follows: 


(Opu/Ox) + (Opv/Oyv) = O 
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the 
components of the mixture, respectively 


where “ and v are v-wise and y-wise velocity 


If we define a stream function W such that 


Oy /Ov pu, Ov /Oxr — pi q 


the equation of continuity |[Eq. (S)| is automatically 


We shall show later that it is convenient to 


use the stream function y 


satisfied. 
as an independent variable 
instead of y which is taken as the distance perpendicu 
lar to the axis of the jet, the v axis 


(c) Equation of Diffusion 


The equation of diffusion is obtained from the con 


servation of mass of one constituent of the mixture 
The diffusion equation consists of three parts the 
diffusion due to convection, the term due to molecular 
transfer, and the term due to chemical reaction. For 


the term due to molecular transfer, we may express it 


in the form of coefficients of diffusion as follows: 


molecular diffusion term 
+ D, grad (In p) - 


Dry grad (In T 


grad [Dp grad « 
(10) 
In Eq. (10), we assume that there is no body force in 
the present problem. 
diffusion, and D, and D, are the coefficients of diffusion 


D is the ordinary coefficient of 


due to gradients of pressure and temperature, respec 


tively. 

Now, for the jet mixing problem, we may apply the 
boundary-layer approximations—i.e., 0/Oy 0/Ox. 
Furthermore, we assuine the pressure is a constant. 
Hence Eq. (10) becomes 
molecular diffusion term 

(0/Oy) (Dp (Oc/Oyv) — (Dr/T) (OT /Oy (10a) 


We assume that chemical reaction takes place ac- 
that the 


cording to first-order chemical kinetics and 


temperature dependence of the reaction rate follows 
the Arrhenius law—1i.e., the rate that the gas | is trans 


formed into gas 2 is 


(l/r 


my Ny 


1/R1 ae 
““-s known as the Arrhenius factor, 7 


the 


energy 


where (1/7)¢ 


time associated with trans- 
the 


energy 


is the characteristic 
and A 
the kinetic 
possess in order to react successfully. 


activation which 


that a 


formation, is 
molecule must 


It is clear that 


Ineasures 


the process could be extended easily to higher order 
single-step reactions; extensions to reactions involving 
several elementary steps is a more serious matter and 
is not covered by the present analysis. 

For two-dimensional steady flow, the diffusion equa- 


tion under the above assumption is 


[O(mycNu)/Ox|] + [O(mecNv) /Oy 
(0/dy) |Dp(dc/dyv)] — (0/d¥) |(Dr/T)(0T/dy) |] — 
(1/r)mNye~ 4/8" (11 
where the terms on the left-hand side of Eq. (11) are 


those due to convection. 
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With the help of the equation of continuity [Eq. 
(8) ], Eq. (11) may be written as follows: 


[m,N/(Be + 1)][u(dc/dx) + v(dc/dy) } 
(0/dv) |Dp(dc/dy) | — (0/dy) [(Dr/T)(oT/dy) |] — 
(1/r)m,Nce~ 4/8" (12) 


It should be noted that the coefficients of diffusion 
and Dr, are, in general, functions of both ¢ and 7. 
In terms of the variables x and y, Eq. (12) becomes 


(dc/Ox) | p?/(N2mymz) |(0/Ow) |p2uD(dc/dw) | — 
[ p?/(N?mymz) |(0/OW) [pu(Dr/T) (oT /dyv) | — 
[c(Bc + 1)/ruje~*/®" (13) 


(d) Equation of Motion 


If the effects on the stress tensor due to diffusion are 
neglected, the equations of motion for the mixture are 
identical to the Navier-Stokes equations. The terms 
of the stress tensor due to diffusion disappear in the 
steady state.’ In the present investigation, we are 
interested in the steady case only; the equation of mo- 
tion will be exactly identical to the Navier-Stokes equa- 
tion. For the two-dimensional steady flow with the 
boundary-layer approximations and constant pressure, 
the equation of motion becomes 


pu(ou/Ox) + pv(du/dv) = (d0/dy)[u(dou/dy)} (14) 


where uy is the coefficient of viscosity which is a function 
of both the temperature 7 and the concentration c of 
the mixture. 

In terms of the variable x and y, Eq. (14) becomes 


Ou /Ox (0/OwW) [upu(ou/dy) | (15) 


(e) Equation of Energy 

Under the same assumptions as those of the equations 
of diffusion and of motion, the equation of energy of 
the mixture in the problem of two-dimensional jet 
mixing of steady flow is 


(0/dyv) [A(O7/dy) ] + 


pu(OC,7/dx) + pv(d0C,7/dy) 
1/R7 (16) 


u(Ou/Ov)? + g(mN,/r)e 


where C, is the specific heat of the mixture at constant 
pressure and X is the coefficient of heat conduction. 
Both C, and \ are functions of temperature and con- 
centration of the mixture; g is the heat release per unit 
mass due to chemical reaction. 

In terms of the variables x and y, Eq. (16) becomes 


(OC,T/dx) = (0/dv) [Apu(oT/ow) | + 
upu(ou/dv)? + g(k/ruje 


, R7 = 
aes (17) 
(III) GENERAL SOLUTION OF Two-DIMENSIONAL 
LAMINAR JET MIXING OF TWO COMPRESSIBLE 
FLUIDS 


From the above analysis, we see that, in the problem 
of two-dimensional laminar jet mixing of two compres- 
sible fluids, there are five unknowns—temperature 7, 
density p, concentration c, and velocity components u 
and v where the pressure p is assumed to be constant in 


this problem. These unknowns could be obtained by 
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solving the five fundamental equations, Eqs. (6), (8 


(13), (15), and (17) with some initial conditions sych 


as: atx = 0 


Cc = oy), T 


u= u(y), v 


Ti(W), p 
vo(W) 


Since we use y as one independent variable, Eq. (8) js 


poly) | , 
¢ 


automatically satisfied. Eq. (6) gives a relation of p 


and 7’in terms of c. In actual analysis, we shall solye 
only Eqs. (13), (15), and (17) simultaneously for the 
unknowns c, u, and 7 with the initial conditions of Eg 
(18). 

One interesting point is that Eqs. (13), (15), and (17 
are the generalized heat-conduction equations. They 
can be solved by a stepwise numerical procedure, finite 
difference method from the initial conditions of Eq. 
(18) without any difficulty because, from the initial 
conditions, it is possible to calculate the rate of change 
of c, u, and 7 with respect to x from Eqs. (13), (15), 
and (17), respectively. By the finite difference method, 
the values of c, u, and 7 at x = Ax may be obtained. 
Hence, by stepwise numerical procedure, the whole 
flow field can be obtained. One advantage of this 
method is that no further assumptions other than those 
already mentioned are required to obtain the solution, 
For instance, we can use any law of variation of vis- 
cosity coefficient with temperature without introducing 
any difficulty in the solution. If one wants to solve 
the problem in a more general manner the computing 
work, of course, is considerable. However, if the 
modern high-speed machine is used for this analysis, 
the work will be of reasonable amount for such an 
analysis to be performed. 

After c(x, ¥), T(x, ¥), and u(x, w) are obtained, the 
y coordinate for any given x may be obtained by the 


following formula: 


(19 


ae 


o pu 


« 


where \ is a constant corresponding to y = 0 at a given 


b 
(IV) SIMPLIFIED SOLUTIONS 


For a general consideration, the problem is ver 
complicated because the density p, coefficient of vis 
cosity u, etc., are not only functions of temperature / 
but are also functions of concentration c. However 
the variations with concentration in many practical 
problems are small, and its effects may be taken int 
consideration by certain mean values. Hence, in order 
to bring out the essential feature of the present prob 
lem, the following simplified assumptions are made: 

(a) The molecular weights of the two components 
of the mixture are assumed to be the same and equal t 
the average actual molecular weights of these two gas¢s 
Then, 

m, = M2, = Mm 
B=0 (20 
p= mN 
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LAMINAR JET MIXING OF 


Hence, the density of the mixture is independent of 


the concentration under this assumption. The equa- 
tion of state, Eq. (6), may be written as follows: 

Pp = p/p T/T 1/T (21) 
where p and 7 are the nondimensional density and 
temperature of the mixture, respectively. Subscript 0 
refers to the values at certain reference conditions. 

b) The coefficient of viscosity is proportional to the 


absolute temperature —1.e., 
a = w/w B(T/T,) = BT (22) 


where B is a constant so chosen that the relation, Eq. 
22), will be a good approximation to the actual law of 


variation of » with respect to 7. Then we have 
pa B = constant (23) 


(c) The specific heat at constant pressure is a con- 
stant—1.e., 
C, = C,, = constant (24) 


(d) The Prandtl Number P, is a constant: 
P, = C,u/X = constant (25) 
(e) The Schmidt Number 5, is a constant: 


S. = p/pD = constant (26) 


c 


(f) The thermal diffusion is negligible: 
D,y = 0 (27) 
We further introduce the following nondimensional 
quantities: 
G=u/U, 2=x/L, ¥ = ¥/WLUBpw = | 
(y/LUm) V R./B, Re = LUpo/mo, Ky (28) 
L/rU, T, = A/RT, K,=4/Crle Ku= 
U2/C,T) = (y — 1)M2, M = U/V yRT)/m 


where L’ and LZ are the characteristic velocity and 
length, respectively, 7 is the Mach Number of the 
flow corresponding to the velocity Ul’, and R, is the 
Reynolds Number of the flow. 

With the above assumptions and the definitions of 
the nondimensional quantities |[Eq. (28)], Eqs. (13), 
(15), and (17) become, respectively, 

d/dx = (1/S.)(0/dv) [a(dc/oP)] — Ky(c/a)e—Te/T 


(29) 
07/0F = (0/dW)[a(Ou/dy) | (30) 


o7/dz = (1/P,)(0/dp) [H(OT/ayP) ] + : 
Kyii(du/dv)? + K,K,(c/a)e—Te/T (31) 


One interesting result is that Eq. (30) is of the same 
form as that of jet mixing of one incompressible fluid. 
Hence, in terms of ¢ and y, the velocity distribution is 
the same as that of an incompressible fluid. However, 
since the temperature and the density distributions are 
different from those of the incompressible fluid flow 
case, the velocity distributions in x and y coordinates 
will be, in general, different from those of the incom- 
Pressible fluid flow case. 
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The present problem depends on the following seven 
nondimensional parameters: 

(a) Reynolds Number R,. This parameter is in 
cluded in the definition of Y. Eq. (19) may be written 


as 


y Yo _ R, [ dy 
i £2 ea 


Q pu 


(b) Prandtl Number P,. 

(c) Schmidt Number 5S,. 

(d) Mach Number .V/ or its corresponding param 
eter Ay. If Ay, is small, the viscous dissipation in 
Eq. (31) may be neglected. 

(e) Heat-release coefficient A, which determines the 
rate of heat release due to chemical reaction. 

(f) Reaction-distance coefficient A, which determines 
the relative reaction distance in the system considered. 

(g) Reaction temperature T,. 

After these seven parameters are given, we may solve 
Eqs. (29)—(31) for #, c, and T from the given initial 
conditions. In general, we may solve Eq. (30) first 
from the initial conditions. Then we solve Eqs. (29) 
and (31) for c and T from the known solution # (#, J) 
and the initial conditions. 

In practical work, one usually wants to solve some 
special cases of the jet mixing problem. A few simple 
cases are given below: 

(1) Jet Mixing of One Compressible Fluid—tIn this 
case, the concentration c is constant. Without loss of 


generality, we may put c = k = 0. We need to solve 
Eqs. (15) and (17) or (30) and (31) simultaneously 
with the initial conditions of Eq. (18). The problem 


has been studied by the author.'. One important result 
is that the properties of the jet mixing depend mainly 
on the momentum of the jet. The effects of the change 
of momentum due to the change in velocity and the 
change in temperature—i.e. density——are the same. 
For a cold jet, the axial velocity decreases more slowly 
than that for a hot jet. The spread of jet is wider for a 
hot jet than for a cold one. Fig. 3 shows three typical 
velocity profiles of two-dimensional jet mixing at dif- 
ferent temperatures ratio of the jet 7 at the exit of the 
nozzle to that of the surrounding stream 7). Ll is the 
velocity difference in the jet from that of the surround- 
ing stream, and Lj) is that value of U, at the exit of the 
nozzle. 

(2) Isothermal Jet Mixing of Two Compressible Fluids 

This is the case of mixing of two gases at low Mach 
Number without chemical reaction. We need to solve 
Eqs. (13) and (15) or (29) and (30) simultaneously by 
neglecting the term of chemical reaction in Eq. (13) or 
(29). This case has also been studied by the author.‘ 
The typical distributions of velocity # and concentra- 
tion c are shown in Fig. 4. The spread of concentration 
is wider than that of velocity. 

(3) Isovel Jet Mixing of Two Compressible Fluids—In 
this case, we consider two streams of different gases of 
the same velocity but of different temperature mixing 
together. Here we solve Eqs. (13) and (17) or (29) 
and (30) simultaneously with “# = 1 and the proper 
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initial conditions. This is the case of laminar jet 
combustion. We shall discuss it further in the next 


two sections. 


(V) IsoveL JET MIXING OF Two COMPRESSIBLE FLUIDS 


We shall consider only the simplified equations, 
Eqs. (29) and (31), which will give the essential feature 
of this problem. Putting “7 = | into Eqs. (29) and 
(31), the fundamental equations for isovel jet mixing 


of two compressible fluids are as follows: 


0c/0% = (1/5S,)(0°c/Ov?) — Kyce P./1 (33) 
oT /oe = (1/P,)(0°T/dP2) + K,Kywe-Te/T (34) 


where we neglect the viscous dissipation terms in Eq. 
(34) by restricting ourselves to the consideration of 
flow with low Mach Number—1.e., J < 1. 

If there is no chemical reaction, both Eqs. (45) and 
(34) reduce to simple heat-conduction equations. Then 
the solution c and 7° may be expressed by simple error 
functions [cf. Eqs. (44) and (46)]. The chemical reac- 
tion terms introduce the deviation of c and T from the 
simple error function. In order to illustrate this 
point, we shall consider two simple examples: 

(1) Mixing of a Two-Dimensional Jet from a Nozzle 


The initial conditions for this case are: at # = O 
Y< = constant = 1,c=1, T=1 } 
y>-—-h = —1, ez1,T=1 ia 
- > = = On) 
) > Wo, c= 0, i i —_ Fs ~ ‘ 

oa * c=0,T7=T7,>1 


where 2 is the total mass flux from the nozzle. For 
simplicity, we put J) = 1. Ti isaconstant. We shall 
take 7; = 3.5 in our numerical example. 

(2) Mixing of Two Uniform Streams*—The initial 
conditions for this case are: at ¢ = 0 


oe ( 
T, > 11 


y> 0, c= 1, 
y < 0, c = 0. 


(36) 


To solve Eqs. (33) and (34) for the initial conditions 
of Eq. (35) or (36), we still require numerical procedure 
for large . However, for small ¢, analytic solutions 
may be found. 

(VI) SOLUTIONS IN THE INITIAL STAGES OF JET MIXING 

Now we introduce a new variable 

f= 9/vVi (37) 
to replace y. 

Eqs. (33) and (34) in terms of # and ¢ become, re- 

spectively, 


(1/.S.)(07¢/O¢?) + (€/2)(d0c/Of) — #(Oc/Ox) = 


@K,ce—Ts/T (38) 
(1/P,)(0°T/os?) + (¢/2)(0T/0f) — #(OT/dz) = 
—¢ K,K,ce~!o/! (39) 


For small values of ¢ we may express the solutions of c 


* This is the case worked out by Marble and Adamson.‘ 
They considered only the simplified equations in this paper. 
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and T in power series of 7 in the following manner 
c=c'(¢) + 2z¢°(¢) + Tic (Cf) + (40a) 


TOS) + ETO) +... + PTO +... 00) 


= 
II 


Substituting Eqs. (40a) and (40b) into Eqs. (38) 
and (39) and collecting terms of the same power of 7 
we have the equations for c'”’(¢) and T(¢). 


For the zeroth power of 7, we have 


(1 S.)(ac" dt? “— 2)(dc d¢) = (dla 


(1/P,)(2T  /de?) + (¢/2)(dT /de) = 0 (41b 


For the first power of %, we have 


(1/.S.) (dc d¢*) + (¢/2 (dco dc) — cP) = 


- { 
K,c*’e 


(1/P (aq de?) + (¢ 2) (dT » df) — e 
Ta 


—K,K,c\’e 
In general, for the mth power of 7, we have 


(1/.S.)(d?c\”’ d¢*) + (¢ 2)\(dc\" /de) — nc” = 
F.(¢) (48a) 


(1/P,)(@T\” de") + (¢ 2)\(dT'\” /dt) — nT" 
= Ke F C) (48b 


where F,,(¢) are known functions of ¢ depending on 
c's and T's where k < n — 1. 

Since Eqs. (41), (42), and (43) are linear, the solu- 
tions for given initial conditions may be easily found 
For instance, for the initial conditions of Eq. (35), the 
zeroth order solutions are 


is = (1 2) ferf[(V/'S, 2)(fo — ¢)] - 
erf[(V/S./2) (fo +6)]; (44a 
T™ = 7, + (C1 — 1))/2]ferf[(V/P,/2) X 
(fo — O)] + erf[(VP, 2)(Go + OF (44d 
where fo = b/Vi 
2 °3 
erf[y] = e *dé 15 
V7 Jo 


The corresponding solutions for the initial conditions 
of Eq. (36) are 
. (1 9)}] | ow 5 9\-1t (46 
( _ 2), 1 + erfl(V S./2)5]5 s 
T =7%4+ ([(1—-—7T))/2]}}l t+ erfl(VP./2¢); (46d 


For high-order functions c'”’ and T‘”’, the boundary 


conditions are: at¢ = +« 
c™(0o) = T'"(o) = 0 (47 


a= fT, 
The solution for c\(¢) or 7“ (¢) is the solution of the 
following typical equation: 


(d?Y/dn?) + (n/2)(dY/dn) — Y = Gi(n) (48 


® or T?) 9 = VS. ¢ or VP, -¢, and G = 


where Y = ¢ 
PF, or —K,F. 


If Y, and Y2 are two linearly independent solution: 
of the homogeneous equation of Eq. (48), the complete 
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Ure 







































































LAMINAR JET MIXING OF TWO COMPRESSIBLE FLUIDS LO17 
ler = | 
| 
(40a) *} | 1 SI 
~~ ee 4 ee eet 
, J > i io ee | > : 
(40b) 004 a | | Ss Lf 
‘ f ' 
/S. (38) i? ‘ 0 )/ 
er of F > 2 | \ ahaa j 
t a | 
(4la 
Fic. 1. Interferogram of a two-dimensional supersonic jet * <a ie eee “ - 
) (41b pproximately at full expansion. pi/p2 = 0.998, MZ = 1.62 
» = pressure at the exit of the nozzle, p2 = pressure of sur 
rounding medium, and J = Mach Number at the exit of the Fic. 5. Temperature distributions in isovel mixing of a two 
zz dimensional jet 
(42a) Y o 
(42b) Jo, Pa, Mm, | GAS | 08 + 4 } = oe See ee 
Tp etc | =a 
x=06 
ge = 6+ 4 4 Shi ee tot 
_----—— RT upetc — 
(43a) .--—. aa 
pa IXI GION at i } 4s 4 i 
Us, Fe, Ma, alias : i} ] 
= GAS 2 ~ 
. u- T etc 
(43b = — | cia | Vy 
: Fic. 2. Jet mixing of two gases 
ling on 
a | 4 oe i ow 
1.0 | j c 4 ; 3 < 
esolu. § .(| | | i [ = Cold Jet > >= v 
found. ° x No Heat Transfer Fic. 6. Concentration distributions in isovel mixing of a two 
5). the — 1 dimensional jet 
—-—Hot vet 
U, 6 Xl 
Ure 
(44a ; . 
4 
(4b iN 
2 
0 T i. catia miniiin 
(45 -6 oS 6 
Y 
ay Fic. 3. Two-dimensional jet in a uniform stream . , “s 
ditions r 
1.0 Fic. 7. Temperature distributions in isovel mixing of two uni 
form streams 
(46a 
(46b .8- L " 
| 4 - 
indary v VAS 
6 + 4 ae a : 
6 / (ff 
(41 C 6 | 4 4 + ttt + + 
of the ro eee 0 | | i | i 
= | 
(48 ad ie T ] 
1G, = | 0 ; P a 
' oF 3 -2 : O 2 3 4 
-6 -4 =2 0 2 a 6 - 


utions Y : , bos aa 
Fic. 8. Concentration distributions in isovel mixing of two 


nplete ‘IG. 4. Two-dimensional jet of CO in a uniform stream of air uniform streams 








1018 JOURNAL OF THE 


solution® of Eq. (48) is 


>” ViG; *” VoGy 
y = J, dn — Vs | — dn + 
2 | i n 1 i n 
Cy Y; + C2 V2 ( $9) 
where 
W = Yi(dY2/dn) — Y2(d Y,/dn) (50) 
and C,; and C; are arbitrary constants. 
The solutions Y; and Y. may be written in the fol- 
lowing forms: 


V¥, = new * + 2y/a[(1/2) + (n?/4) Jferf(n/2) — 1) 
Yo = ne” 1 + 2 w[(1/2) + (2/4) ][erf(n/2) + 1]8 


(51) 
Then 
W = -¥V re" (52). 
The boundary conditions that at 7 = +o, Y = 0 
give that 
= 4=0 (53) 


and that the limits of the integrals of Eq. (49) should 
be taken in the following manner: 


Y2G, VG, 
€ 
« v4 u] Yr 


The expression, Eq. (54), may be evaluated by numeri- 


y= Y e”/4dn + Vo Sdn (54) 


cal integration. 


(VII) NuMERICAL EXAMPLES OF ISOVEL JET MIXING 


Numerical examples of isovel jet mixing of two com- 
pressible fluids have been worked out based on the 
following data :* 


Fr, = O91, - = 1.00, K, = 4.497 t (55) 
K,=5xX 10", # = 1.00, T, = 83.870f ‘°° 
The initial conditions are: 
(a) Two-Dimensional Jet—At = = 0 
y<landy> -l: c=1, T= | (56) 
a a 06 
y>landy< —-l: c=0, fT = 3.5) 
(b) Mixing of Two Uniform Streamst—At % = 0 
y>0: c=1, T= 1.00\ (a2) 
y<0: c=0, T=3.50S - 


Eqs. (33) and (34) have been integrated numerically 
based on the values of the parameters of Eq. (55) and 
the initial conditions of Eqs. (56) and (57) by the use 


* These data are the same as those used by Marble and Adam- 
son.* 


t This case has been solved approximately in reference 3. 
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of a high-speed calculation machine.{ The computa. 
tion is carried up to a point where the essential features 
It should be 
noticed that the concentration ¢ is a nonnegatiy 


of the solution may be clearly indicated. 


quantity. 
ped out whenever it has the tendency to give negatiy, 


The chemical reaction term should be drop 


value of c. 

Figs. 5 and 6 give the temperature and concentration 
distributions, respectively, for the initial conditions of 
Eq. (56). Fig. 5 shows that there is a tendency to «& 
velop a flame front at the edge of the jet—i.e, 
peak in temperature distribution. 

Figs. 7 and 8 give the temperature and concentration 
distributions, respectively, for the initial conditions oj 
Eq. (57). Fig. 7 also shows that there is a tendency t 
develop a flame front at the boundary of mixing 
These results are qualitatively the same as those given 
by Marble and Adamson. 

Calculations of first-order solutions c”’ and 7" for 
the initial conditions of Eqs. (56) and (57) have also 
been carried out. Qualitatively, the first two terms oj 
Eqs. (40a) and (40b) will give the same distributions 
as those obtained by numerical integration, but, quanti- 
tatively, the results of the first two terms give tempera- 
ture distributions higher than those obtained by the 
numerical integration. 

The coordinate y in Figs. 5—-S may be easily changed 
into y by the help of Eq. (32). 
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Thermal Buckling of Supersonic Wing Panels’ 


N. J. HOFFi 
Polytechnic Institute of Brooklyn 


SUMMARY 


The temperature and thermal stress distributions are analyzed 
in multicellular supersonic wing structures. A buckling criterion 
is established for the panels of cover plates subjected to thermal 


stresses 


INTRODUCTION AND SUMMARY OF RESULTS 


_— THE SURFACE of a multicellular supersonic 
wing is heated by the air flowing in the boundary 
layer (see Fig. 1), the temperature of the cover plates 
rises and the cover plates expand. This expansion is 
resisted by the shear webs which remain comparatively 
cool as they are not in contact with the boundary layer 
and are heated only by the cover plates through con- 
duction and possibly through convection and radiation. 
The resistance of the shear webs sets up compressive 
stresses in the cover plates and tensile stresses in the 
shear webs. Under the compressive stresses so arising, 
the cover plates have been observed to buckle. 

Buckling of the cover plates is objectionable not only 
for structural reasons but also because the buckled, 
wavy surfaces alter the flow of air, lead to losses in 
aerodynamic efficiency, and may even cause aerody- 
namic disturbances resulting in failure. It is impor- 
tant, therefore, for the airplane designer to be able to 
predict the onset of thermal buckling—that is, buckling 
caused by stresses arising from a nonuniform heating 
of the structure. In multicellular wing structures, the 
variation of the temperature through the thickness of 
the cover plates is usually not pronounced, at least 
not when the speed of flight does not exceed Mach 
Numbers of 3 or 4. Similarly, with a wing of normal 
aspect ratio, the variation in the spanwise direction 
may be disregarded in order to simplify the analysis. 
However, the heating rate generally varies noticeably 
in the chordwise direction and with time. 

As all possible variations in flight path and construc- 
tion cannot be considered simultaneously when formulas 
to be used in practical work are developed, it is assumed 
in the present analysis that, over one cell of the wing 
and during the period of time considered, the heat in- 
put g from the air to the cover plates is constant. In 
in actual problem when gq is a variable, substitution of 
its average value into the formulas here derived should 
give an indication of the danger of buckling. Heat 
transfer from cover plate to web through the interior of 
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the wing by means of radiation and convection is dis- 
regarded; it is often unimportant as was discussed by 
the writer.'~* The only means of heating the web is 
therefore conduction, and the assumption is made that 
no resistance to the flow of heat is offered by the con- 
nection between cover plate and web. When poor 
thermal connection diminishes this heat flow, an ef- 
fective thickness h,, less than the geometric thickness 
may be assumed for the shear web in the analysis 

In the calculation of the thermal stresses, the shear 
webs are assumed to be flat plates rigidly connected to 
the cover plates along their lines of intersection. The 
thermal stresses can be diminished by providing non- 
rigid connections. The final formulas of this paper 
are, of course, not valid for wings incorporating such 
arrangements. 

In the analysis of buckling, the cover plates are as- 
sumed to be perfectly flat, to have a width L,, and to 
be long in the spanwise direction. Without this last 
assumption, it is not justifiable to consider the thermal 
stresses constant in the spanwise direction. Of course, 
they vary in the chordwise direction. It is further as- 
sumed that conditions change sufficiently slowly from 
cell to cell so that two or more adjacent cells reach their 
stability limit at about the same time. This means 
that no portion of the cover plate is noticeably re- 
strained in its buckling by the adjacent portions of the 
cover plate. If the connection between cover plate and 
web is flexible enough, each panel of the cover plate 
can be considered as a rectangular plate of width L, 
simply supported along the webs. 

It is hoped that the various equations presenting the 
laws of temperature and thermal stress distribution in 
multicellular supersonic wings, as well as the formulas 
from which the buckling load can be calculated for dif- 
ferent thermal stress distributions, will be of some 
practical help to designers of supersonic aircraft. The 
use of these formulas must be adapted to the individual 
design and flight conditions. However, some simple 
results may be of sufficiently general interest to be 
quoted in this introduction. They are rigorously valid 
only under the highly idealized conditions for which 
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they were derived, but they can provide an estimate 


even under conditions differing widely from the ideal 
ones. The first result of this kind concerns the tem- 
perature distribution in the cover plates and the second 
the buckling load. 

If a series of identical cells, such as the one shown 
in Fig. 2, is subjected to aerodynamic heating of a con- 
stant intensity g (B.t.u. per sy.in. per sec.), the teim- 
perature of the cover plates, in the absence of shear 
webs, would rise in accordance with the equation 


T° Ht (1) 


where 7’, is the umform temperature of the cover plates, 


tis time, and 
[1 q h.cpg (2) 


with /, the thickness of the cover plate, c the specific 
heat of its material, and pg the weight density of its 
material. Under most practical conditions, this is also 
the temperature in the cover plate midway between two 
adjacent shear webs in very good approximation. In 
the derivation of this formula it was assumed that the 
two cover plates were heated equally and that no heat 
was lost by the cover plates through radiation or con- 
vection. 

At the juncture with the shear webs, the temperature 
of the cover plates is lower in consequence of the heat 
flowing into the web. There the temperature is given 
by the equation 


Teo = (1 — p*)Ht (3) 
where 


(1/p*) 1+ (Ae/hty) (Ke/ Kw) (Qe/ae) (4) 


K is the conductivity, a is the diffusivity of the mate- 
The subscripts c 
When 
cover plate and web are of the same material, Eq. (4) 


rial, and / is the plate thickness. 
and w refer to cover plate and web, respectively. 


reduces to 


(1/p**) 1 + (h./hy) (4a) 


The second interesting result refers to plates in which 
the thermal stress distribution can be represented 
approximately by a uniform compressive stress upon 
which is superimposed a second harmonic of the com- 
pressive stress whose amplitude is p times the average 
stress. The buckling stress o,, of the long rectangular 
plate of width LZ, and thickness h, is given by 


1d (5) I 
: - (5) 
3(1 — »,”7) \L./ 1 — (p/2) 
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where /, is Young's modulus, y, is Poisson's ratio of th, 
material of the cover plate, and / is the ratio of th 
second harmonic of the compressive stress to the ayer 
this formula, o 


age compressive stress. In is th 


average compressive stress—that is, the total compres 
sive load acting on the plate divided by L,h, eve 
though the load is distributed nonuniformly as ind; 
cated in Fig. 3. In the development of this formulg 
the plate was assumed to be simply supported along 
the shear webs. Eq. (5) is accurate enough for pra 
tical use when —2 < p< 1. 

Finally, it should be noted that a plate cannot by 
expected to remain perfectly flat until the critical load 
is reached and then to develop large deformations sud 
denly. Because of the unavoidable deviations from 
perfect flatness and in consequence of variations in 
thickness and material properties, a real plate, as dis 
tinct from an ideal one, begins to deform as soon as 
thermal stresses develop in it in consequence of the 


heating. When the theoretical critical condition js 
reached, the deformations increase rapidly. This jis 
the physical significance of the buckling load. The 


buckling condition can be determined in an experiment 
if the test results are analyzed by means of a Southwell 
plot or if such criteria as the top-of-the-knee and the 
inflexion point are used. 

As far as is known to the author, the only earlier pub- 
lication dealing with thermal buckling was written by 
Gossard, Seide, and Roberts‘ in 1952. In this paper the 
postbuckling behavior of a rectangular plate of finite 
side ratio was analyzed; the temperature distribution 
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THERMAL BUCKLING OF 


«as assumed to be triangular in the chordwise and 


onstant in the spanwise directions.f 
CALCULATION OF THE THERMAL STRESSES 


Temperature Distribution 


The wing structure is assumed to consist of two 
parallel plane cover plates connected by a series of 
equidistant webs (see Fig. 2). The length of the 
structure perpendicular to the plane of the figure is 
The thermal connection between any two 
elements is perfect. The cover plates are uniformly 
heated at the rate of g B.t.u. per sq.in. per sec. If the 
webs were not heated by conduction and if no heat 


verv great. 


were lost by convection and radiation, the temperature 


{f the cover plate would increase at the rate 


oT /dt = g/heepg = H 6) 


where ¢ is time, 7’ the temperature of the cover plate, 
the specific heat of its material, p.g the weight 
lensity, and /, the thickness of the cover plate. Of 
course, heat is lost to the web by conduction. For 
this reason, the temperature is nonuniform over the 
structure, and the temperature distribution in the cover 
plate is governed by the differential equation 
OT /dt = a,.(0°T/Ox") + (¢/hCcpc£) 7) 
where a, is the diffusivity of the material and x is the 
distance measured along the cover plate from the junc- 
ture between cover plate and web as shown in Fig. 2. 
As the web is heated only by conduction, its tem- 
perature must be calculated from the equation 


a,,(0°T'/O& 


(S 


oT /ot 


where € is the distance measured along the web from 
the juncture between cover plate and web as shown 
in Fig. 2. Eqs. 7 and 8 are well known; they can be 
found, for instance, in reference 5. 

At the juncture of cover plate and web, the temper- 
iture of the cover plate must be equal to that of the 
web. Also, the heat lost by the cover plate must be 
equal to the heat gained by the web. Hence two 
boundary conditions can be given as 

i T 


when x =& 0 (9) 


P, _ ol, — P; - oh,/2 (10) 


x 


itall times. In the last equation, ® is the flux of heat. 
Itis given by 
@ = —K(dT7/Odx) (lla) 


in the cover plate and by 


? —K(OT/0&) (11b) 


in the web. 
the material. 


t After completion of this manuscript, another paper dealing 
vith the subject was published by Joseph N. Kotanchik, Aldie 
E. Johnson, Jr , and Robert D. Ross under the title Rapid Radi- 
int Heating Tests of Multiweb Beams, NACA TN 3474, September, 


The letter A denotes the conductivity of 
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The solution of the temperature distribution prob 
lem is made easier if one cover plate and one web are 
assumed to extend from their juncture to infinity and 
if the second cover plate and all the other webs are 
disregarded. It will be shown later that this approxi 
mation is a satisfactory one under many practical con- 
ditions. If this assumption is made, two additional 
boundary conditions are that, far away from the junc- 
ture, the temperature of the web must remain un- 
changed and that that of the cover plate must rise as 
if no webs existed. Finally, one may assume that 
at ¢ 0, when the heating begins, the temperature of 
the structure is uniformly zero. 

A particular solution of Eq. (7) is 


z- Ht 12 


as can be easily proved by substitution. This solution 
represents the temperature history of the cover plate 
in the absence of webs. 

It can now be conjectured that the temperature of 
the web-cover-plate juncture will also rise linearly with 
time—although with a proportionality factor p rather 
than //. If one end of a fully insulated semi-infinite 
rod is heated in such a manner that the temperature 
of the end x O is pt, the temperature distribution in 
the rod is, according to page 45 of reference 5, 


r tpt rerte(x/2 Yat 13 


where 7*erfcz is the twice iterated complementary error 
function of z discussed in some detail in Appendix 2 of 
the same reference. As, in the present problem, heat 
is taken away from the cover plate at x = 0, p must be 
negative with a yet unknown value. The complete 
solution can be written as 

T(x, t) tpt teric(x/2~7/at) + Ht 14 
This solution obviously satisfies the boundary condi 


tion atinfinity. At the juncture it yields 


T (0, t) p+H)t 15 
Differentiation gives 
O7'/Ox —2p(t/a,)" “ierfe(x/2VY acd 16 
which becomes at + = 0 
(O7/Ox), —2p(t/ma,)' * 17 


Hence, at the juncture, the flux is 


®,_, = —K,(OT/Ox), 2K .p(t/xa,)'* IS 

If the rate of change of the temperature in the web 
at £ = Ois7, the temperature in the web is 

T(é, ¢ trt rerfe(E/2Y dy 19 


The temperature and the flux at the juncture are 
T.,,(O, ¢ rt 20 
(0, ¢ 2K_r(t/ra,)'* 21 


The requirement of equal temperatures at the junc- 


ture can be satisfied only if 
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(22) 


p+H=r 


As each web receives heat from two cover plates, 
the equality of heat lost and that gained can be ex- 


pressed as 


2K .p(t/ma.)"*h —2K yr (t/may)*(Ny/2) (23) 


From these equations, » can be calculated as 


p= —ptH = —H/[l + 2(h-K/ipKy) (dy/ae)*] (24) 


When the material of cover plate and web is the same, 
p* becomes p** defined as 


(25) 


1/p** = 1 + 2(h./Nw) 


With this expression given, the temperature distribu- 


tion in the entire structure is known. It must be 
determined, however, under what condition this solu- 
tion is acceptable in view of the simplifying assump- 


tions. 


Range of Validity of Assumptions 

The assumption of an infinite length for cover plate 
and web in the x and & directions means physically that 
the temperature distribution in the neighborhood of 
a cover-plate web juncture is analyzed without regard 
to the effect of the junctures of the other webs with the 
cover plate and of the same web with the second cover 
plate. In an engineering calculation, such a procedure 
is certainly permissible if it can be shown that the tem- 
perature in the cover plate midway between adjacent 
webs and that of the web midway between the two cover 
plates is essentially unaffected by the heat flowing 
through the junctures. This can be stated more ex- 
plicitly by stipulating that the value of 7’erfe(v/2Y_ .t) 
in Eq. (14) and that of /erfe(é 2+/a,t) in Eq. (19) 
should not be greater at the middle of the panels than 
about 5 per cent of the value at the juncture. 

As 47°erfc(0) = 1 and 47°erfe(1) = 0.0568, the re- 
quirement stipulated is satisfied if, in the middle of the 
panels, the argument of the error function is unity or 
greater. In the form of an inequality 


(26) 


L ‘44/at ee 


(27) 


or L?/16at > 1 


When the diffusivity of the material and the length of 
the material and the length of the panel are given, 
this is a restriction on the time during which the 


formulas are valid. It must be required that 


t ee L?/16a (28) 

A few typical values of the diffusivity are: 
0.0075 in.*/sec. for 18-8 stainless steel 
0.0186 in.?/sec. 
U: 130 


for structural steel 


in.*/sec. for pure aluminum 


If the webs are spaced 10 in. apart, the inequality, 
Eq. (28), is satisfied for the cover plate if 


t 833 sec. for 18-8 stainless steel 


< 
t < 336 sec. 
< 46° sec. 


for structural steel 
for pure aluminum 
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These time limits are broad enough to include gen 
erally the most important period from the standpoint 
If it js 
necessary to extend the period of time beyond th 


of thermal stresses and thermal buckling. 


limits given, the basic equations are still valid, but th, 
boundary conditions at infinity must be replaced by 
boundary conditions at a finite distance. The tech 
niques to be used in such calculations are discussed jy 


references 5 and 6. 


Thermal Stresses in an Infinite Plate 


The nonuniform temperature causes a nonuniform 
expansion of the material. The nonuniform expansion 
following directly from the law of thermal expansion 
would be possible only if cracks developed in the 


structure. As no cracks develop—at least not until 


the temperature differences become large—they must 
be prevented from occurring by deformations caused 
by internal stresses. In a long cover plate of width 
L, these internal, or thermal, stresses can be calculated 
on the basis of the following considerations: 

Let us assume that the plate is cut up into a large 
number of narrow strips in the longitudinal direction 
before heating begins. Next, heat is applied to the 
strips in such a manner as to obtain in them the tem- 
perature distribution calculated earlier. As the middle 
of the plate is the warmest, it expands more in the 
longitudinal direction than the rest of the plate. The 
narrow edges of the long plate do not therefore remain 


straight. The elongation of each strip is 


AL al(x, i)L" 29 


if a is the coefficient of thermal expansion and L* 
is the length of the plate. The average elongation of 
the panel is 
*/ 
AL... (ai /L) T(x, t)dx 3U 
Let us now apply to the ends of the strips stresses of 
the magnitude 
¢ = ak(T — T,,, 31 
with the positive sign indicating compression; during 
the application of these stresses the strips still remain 
disconnected. The stresses cause strains of a magni- 
tude 
€ a(T — T, be 
where again the positive sign indicates compression 
The shortening 
AL = aL*(T — Tare) 30 
superimposed upon the thermal elongation re-establishes 
the straight-line shape of the shorter edges of the plate, 
but the initial length L* 
amount 


is uniformly increased by the 


AL® = al Tas 34 


in consequence of the average thermal expansion. 
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At this stage of the process, the connection between 
adjacent strips is re-established. By this sequence of 
events, a plate is obtained in which the internal strain 
is uniformly zero, the temperature varies across the 
width, and internal stresses exist in accordance with 
Eq. (31 
the width, but they are constant in the longitudinal 


The internal, or thermal, stresses vary across 


direction; at the shorter edges they are in equilibrium 


with the stresses applied there. 


It is easy to see from the formulas that these applied 
stresses add up to zero resultant. For this reason, they 
can be removed without influencing the internal stress 
distribution anywhere in the plate except in the end 
zones extending to a distance of about LZ from the free 
edges (if L is the width of the panel). 
quence of St. Venant’s principle. But, with the ex- 
ternally applied stresses removed, the plate is in equilib- 


This is a conse- 


rium under the thermal stresses just calculated, and 
its deformations are compatible without the formation 
of cracks. Hence, in the major portion of the plate, 
the thermal stress distribution is as given by Eq. (31). 


When the temperature of the plate is not symmetric 
with respect to x = L/2, a term of E7*x must be 
added to the right-hand member of Eq. (31). 


Thermal Stresses in the Wing Structure 


These considerations need be modified only in minor 
detail if the thermal stresses are sought in the super- 
sonic wing structure. Because of the double sym- 
metry of the idealized structure of Fig. 2, it suffices to 
consider one half of a cover plate and one half of a web, 
with the idealized web having only one half the thick- 
ness hi, of the actual web (as two cover plates feed 
heat into one web). This is true provided the two 
cover plates are heated in the same manner. 
assumed that the inequality, Eq. (28), is satisfied. 
This means that the temperature distribution derived 
for the semi-infinite rods is a satisfactory approximation 
to the actual temperature distribution in the structure. 


It is also 


Under these conditions, the average temperature of 


the cover plate is 


) Ple/2 
1 — T(x, t)dx | 
Lose as 
(35) 
9 ¢ 9 fs 
T(x, t)dx — | T(x, t)dx 
£ J Lew left 
¥ ; ; a 
As al’eric Bx dx 3 Merfc Bb 
Jb i. 
and since 7*erfc(0) 1/(OV mw) 0.0940 
one obtains 
Tenis Ht \1 — (2p*/L,) (at)? [0.752 — 
Si®erfe(L/4 ~/a,t)|{ (36) 


In many practical applications the last term is 


insignificant. When L/4\/a, = 1, the quantity 
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Hence, omission of the last 
For this 


Sieric(L/4+/a) = 0.0291. 
term results in an error of only 4 per cent. 
reason, it is recommended that the average temperature 
of the cover plate be computed from the formula 


ee HIt }1 — 1.5(p*/L,) (ad)"*} (37) 


This formula is naturally valid only if the inequality, 
Eq. (28), holds. 
The average temperature of the web is 


: eee? 1.5Ht((1 — p*)/Ly| (Aut) * (38) 

The thermal stress in cover plate and web is 
G, a,EHt 1 — 4p*iPerfe(x/2 ~V/a,t) — M; (39) 
On a, dit }4(1 — p*) Perfe(E/2 Yayt) — Mj (40) 


where the constant ./ must be determined from the 
equilibrium condition 


L-/2 eLw/2 ] 
1 
f oh. dx + } Ow — dé 0 


(41) 


| fas 
(42) 
3 T 


and with a similar expression for the web, integration 
and evaluation of .V/ yield 
(1 + B)M = (8/3) (1 — p*) (ayt/aL,?)"* + 


Bil — (8 3)p* (act rL,*)' "7 (43) 


where B = 2a Ae) eu he) | (44) 
or approximately 
M = [1/(1 + B)] $1.50. — p*) (dut/Ly?)"* 4 

Bll — 1.5p*(a.t/L.2)"*}{ (45) 


When the material of cover plate and web is the same, 


M is replaced by ./* defined as 


Ae Pes & (46) 


(h,/h-) 


Representation of Thermal Stresses by Fourier Series 


In the evaluation of the danger of thermal buckling, 
the thermal stresses acting in the cover plate must be 
expanded into a Fourier cosine series. The first term 
of the expansion is the average stress which can be 


easily found as 


Saree = Geli Ht {1 — M — 1.5p*(a.t/L,2)"*} (47) 


In the calculation of the higher harmonics, the following 
integral must be solved :f 
Carrier for the 


+t The author is indebted to Prof. George F 


evaluation of this definite integral. 
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Fic. 4. Graph of function (A) in Eq. (59) 


pore NWN - x 
Sy = cos tt 7°erfe == 0% 
v0 ZL, 20V/ at 
( {S\) 


{ MEX yt, —* _ g 
— cos “eric ax 
0 iE, 2V/at 


This infinite integral will be evaluated with the aid 
of the Laplace transformation method. The Laplace 
transform @(p) of a function v(t) is defined as 


a(p) | e~”'v(t) dt (49) 


According to page 3S0 of reference 5, the Laplace trans- 


form of the function 


o(t) = 4t Perfe(x/2V/at) (50) 
is a(p) = (1/p?)e-*(e/a)/* (51) 
Consequently, 
- n7x | 1/2 : 
J, = vind Lp? e—x(b/a)"’* dx (52) 
0 4 - 
As cos (nmx/L) = Re e'"**/* (53) 


one may write 


l 9 
- re f  ex—(/a)'/2 + in (e/L)] as| 
0 p- 


ex[—(p a)l/2 + i(nx/L)] } | 


I 
= Re 4— mediates , 
lp? —(p/a)"" + Une L)$, 
At the upper limit the numerator vanishes, while at 
the lower limit its value is unity. Hence 
(1/p*) (1/p?) (p/a)'”” 
—(p/a)’? + i(nn/L)  (p/a) + (nx/L)? 


Jn = —Re 


a? 
p*[p + a(nr/L)?} 


(55) 


This expression can be rewritten as the partial fraction 


.. ._ 3 ~ 1 ( 
n?xa’? p' 2p? Ip + a(nn/L)?]$ 


(56) 


Finally, the inverse transform of this expression jg 
needed. In reference 7 one can find the result 

J, = (2/A2) (at®/r) P(A 37 
where A nat) */L 58 


A 


¢ eA 
and (A) = 1 — } e* dx 59 
A 


The following expansions are useful in computations 
@(A) ~ 1 — (1/2A"*) — (1/4A4) — (3/S8A5) — 
fA>3 (59a 
(A) = (2/3)A2(1 — (2/5)A2 + 
(4/35)A? — (S/315)A®... (59b 
Eq. (59b) should be used when 0 < A < 0.7. The 
function ®(A) is shown in Fig. 4. When A is greater 


than about 5, ®(A) can be taken as unity. The expres 
sion needed in the calculation of the Fourier coefticients 


is 
J, = (0.114/n?)L2(t/a)' *@(A (60 
The nth harmonic 7, of the stress is 
; $  Le/2 MTX 
ty > cos a. ax n 4&0... © 
oe, P L 


Substitution yields 
In = —Q, EH p* J,(4/L, / 
—0.456 a-E.HL,p*(1/n?) (t/a-)' *® (A)\ 
The value is given in psi if the units lb., in., and sec 
are used. 


BUCKLING ANALYSIS 


The Equations Defining the Problem 


The flat plate of uniform thickness / shown in Fig. 4 


is subjected to compressive thermal stresses in the j 
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Fic. 5. Plate subjected to thermal stresses 
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tirection lhe four edges of the plate are simply sup . nT = kr i 
a ' : . at vi = sin : a, sin ov 
norted Che buckling problem is governed by the dif , Lb én 
ferential equation 
Substitution yields the equation 
Vw —(h/D)ow,,, (63 
, : . > oi 2\2 ot kawx A nl > > 
where w is the deflection of the plate perpendicular to LZ, \" k*)* sin G 90 = & 
the plane of the plate, subscripts following a comma 7 
indicate differentiation, the bending rigidity D of the _ (Rk + m)rx _ (kh m) mx | 7 
sin + sin 70 
plate 1S i L f 
D Eh*/{(1201 — v (64 ie nae , 
, Chis condition of equality between the two trigo 
nd the operator V4 is defined as nometric series is satisfied identically only if it is 
; satisfied for each sine term. Consequently, 
Ty V2)? [(0?/Ox?) + (07/dv?) |? w 
2C 2 
O'w/Ox*) + 2(0*w/Ox-Oy") + (O*w/Oy" (65 (n? +7 a ® Pm (Aj4 + §.a (7 
n 
The boundary conditions are 
where the subscript 7 — m is to be taken as an absolute 
P 07w /Ox? Q when vx Oh ¥ ; 
a value 
(Ob 
w = 0°w/Oy? = 0 when y = 0, L/n\ : 
> - a a 
Eqs (63) to (66) are well known; they can be found, for 
instance, in reference S. Moreover, 1 when 7 > m | 
[he compressive stress is represented by the series bin 0 when j m 2 
—1 when 7 mS 
o marx = 
(a - > B Pn cos (04 
Cn - and j = 1, 2, 3, (73 
1 ‘ > ‘ »\ 7t . 
wher a w?E/({12(1 — v*)}i (A/L (OS re ee or rg: ei 
_ t | M Phis infinite set of linear equations in the coefficients 
. a, is homogeneous. It can have a nontrivial solution 
Solution of the Equations ss Bl at na ea , 
q only if its matrix is singular. The condition for this 
rhe solution of the characteristic value problem is is that the determinant formed of the multiplying fac 
issumed in the form tors of the a,, vanish 
2(po — CK,) — p2 i — 2 p2 — Pps 
pi — Pp 2(po — CKe) — ps Pi — Ps 
po — Ps Pi — Ps 2(po — CK3) — p = 
= {) ( }) 
where K, = [(n? + s*)/n}? 
When the p,, are prescribed, Eq. (74) defines the value of the divisor C in Eq. (67) that corresponds to buckling. 
With C known, the ratios 
Qa, = a;y/aQ) (40) 
can be calculated from the set of equations 
(pi — ps)as + (p2 — ps)as + (ps — ps)axs + = —2p, + 2CA, + p 
(2p9 — 2CK2 — ps)a2 + (pi — Psda + (pe Pe) Qs + =— pi tT P 
(pi = Ps) a2 4 (2p, a 9CK ‘oie Pea + (pi — pz) as ale = — pot py - 
(46 


Uniform Stress 


The buckling stress will next be calculated for a few particular cases. 


all the other Pm are assumed to vanish. 
terminant vanishes when 


Then the matrix of Eqs. (71) degenerates into a diagonal matrix. 


As a first example, /» is taken as unity, and 
Its de- 


(1 — CK,) (1 — CK) (1 — CK;) (1 — CKy)... = 0 (77) 
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from which it follows that 


(1/C = K 78 
The smallest buckling load is obtained when 
K, = Ky = [(nm? + 1)/n/? 79 
and the minimum of this corresponds to n = 1: 
K, min = 4 80 


It follows then from Eqs. (67), (68), (78), and (90) that the minimum buckling stress of a rectangular plate subjected 


to uniform compression is 
Ser = 400 = [x°E/3(1 — v?)] (4/L)? 8] 
which is well known. 
Second Harmonic Superimposed on Uniform Stress 
Let us assume next that 
po = 1 be = p Pm = 0 when m = 0, 2 (82 


With these values, Eq. (74) becomes: 


(1 — CK) — (p/2) 0 p/2 0 0 0 0 
0 (1 — CKe) 0 p/2 0 0 0 
p/2 0) (1 — CK3) 0 p/2 0 0 
0 p/2 0 (1 — CK;) 0 p/2 0 
0 0 p/2 0 (1 — CK;) 0 p/2 i ae 
0 0 0 p/2 0 (1 — CKe) 0 ' 
0 0 0 0 p/2 0 (1 — CK;) 


The matrix of the equations can be transformed into a triangular one by the following operations: Multiply the 





seventh row by p/(2(1 — CA;)| and subtract the result from the fifth row. This is permissible if 1 — CA; # 0. 
The operation eliminates p/2 from the fifth row and the seventh column but adds — p?/4(1 — C&A;) to the principal 
diagonal member of the fifth row. Next p/2 is eliminated from the third row and the fifth column by multiplying 
the new fifth row by (p/2)/)(1 — CK;) — [p?/4(1 — CK;)]} and subtracting the product from the third row. This 
is permissible if the denominator is not zero. The new third row is now used to get rid of p/2 in the first row and the 
third column. If a similar sequence of operations is performed on the even-numbered rows, a matrix is obtained in 
which every member to the right and above the principal diagonal is zero. The value of the determinant of such a 
triangular matrix is the product of the principal diagonal members P,. Eq. (83) can therefore be written in the form 
of an infinite product 


YP. = @ (84) 
s=1 
In this equation P, is the infinite continued fraction 
oe (p/2)? 
P, = [1 — CK, — (p/2)] - Pi: - 
. (p ” 
(1 — CK;) — t 
2 (p/2Z)* 
(1 — CK;) — {P 
(l — CK;) —_ Sd 
Similarly, : (p/2) 
° P; = (1 — CK;) - f ste 
: (p/2) 
(1 = CK;) P . 
(1 — CA;) —... (30) 
: (p/2)? 
P, = (1 — CK.) — f 
as (p »\2 
(1 — CK,) sas saa od 
(1 — CK.) —... (d/) 
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Eq. (84) 15, of course, satisfied if any one of the p, 
vanishes. Of practical importance is the highest 
value of C obtainable from this condition because it 
corresponds to the lowest buckling load. 

Let us calculate the buckling load of a square panel; 


then 


AK 25. & 100, A; 289, 
1,369, A; 2.500 


A; 676, Ke 

If P; is to vanish, in a first approximation one may write 

1 — C,A, — (p/2) 1 — 4C, — (p/2 Q (88) 
from which 


C, = (1/4) [1 — (p/2 (S9 


In a second approximation, one may disregard the ex- 
pressions containing A, with s greater than 3. Eq. (S85) 


then reduces to the quadratic 


(1 — 4C, — (p/2)] (A — 100C,) (p/2)° (90) 


whose approximate solution 1s 
C, = C1 + « (91) 
where 
p [(p/2)?]/[24 — 49(p/2) + 25(p/2)?] (91a 
It can be seen that the correction is very small when p 


Even if the amplitude of the second harmonic 
1, the correction 


is small. 
is the same as the average stress, p 
is negligible in engineering calculations—-namely e¢ 
23. Hence Eq. (SS) is accurate enough when —2 < 
p<. 

If P2 or P3 are required to vanish, the first approxi- 
mation yields 0.04 and 0.01 for C; these correspond to 
much higher buckling loads than the ones already ob- 
tained and are unimportant in practical applications. 

On the assumption that the first approximation is 
sufficiently accurate in the range —1 < p < 1, numeri- 
cal values were computed for C and are shown in Fig. 6. 
rhe curve indicates that the total compressive force 
that is, the integral of the stress over the width L of the 
plate times the plate thickness / necessary for buck- 
ling is lowered 33 per cent when the compressive stress 
is concentrated in the middle of the plate in such a 
manner that its value there is twice the average value 
while, at the edge, the stress is zero (p —1). On 
the other hand, when the compressive stress in the 
middle is reduced to zero while that at the edge is 
doubled (p = 1), the average stress necessary to cause 
buckling is doubled. 

The panel can buckle even if the net compressive 
load is zero. If, for instance, the assumption is made 
that 


when s ~ 2 (92) 


transformations of the kind undertaken earlier result in 
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Fic. 6. Variation of C with / 


P —CK, 
~ (1/2) — (1/4 
ac (1/4 
—CK; — 
aie 1/4 
—CK,; — — 
— CK; (93 
sat (1/4 
P» —CK, — 
CK i 
~ i aes 
—CK, —... O4 
The second approximation to the smallest buckling 
load is 
1/C ~ —7.67 (95 
on the basis of ?; and 
1/C 2(K-K,)*"" +170 (96 


on the basis of P». 


Buckled Shape 


In the calculation of the higher harmonics of the de- 
flected shape at buckling, the value of C must be very 
accurately known because small differences of large 
As the 


value was computed for = | from the first two odd- 


numbers occur in the solution of Eqs. (76). 
numbered rows of Eq. (83), it would be inconsistent to 
determine the amplitudes of modes higher than the 
third. If one takes 

by = 1 bo = | C = 3/23 (97 
the first odd-numbered equation of the set Eqs. (76 
yields 


(98 


The panel therefore buckles almost entirely according 
to the simple half-sine-wave pattern even though the 
loading is far from uniform when p = 1. 

When p2 = | is the only nonvanishing p,, buckling 
takes place in a combination of the first and the third 
modes; with the first mode prevailing. This can be 
seen if 1/C = —7.67 is substituted from Eq. (95) into 


the first of Eqs. (76). The result is 


a; = 1/2CK; = —0.0384 (99) 
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Finally, the buckled shape is exactly the half sine 
wave when py) = | and p, = O for all other values of s. 
Then the buckling load is defined by Eq. (Sl). The 
numerical coefficient in this equation was obtained 
from the first of Eqs. (74) alone, and thus no mode 
shape higher than the first can have an influence on 


this type of buckling. 


NUMERICAL EXAMPLE 


To illustrate the application of the theory, a numeri 
cal example is worked out in some detail. The center 
section of the wing is assumed to be represented with 
sufficient accuracy by Fig. 2 with the following numeri- 
cal values: 
= Sin. 
= 0.0625 in. 


b, = 12 om. £ 

h, = 0.125 in. h 
The wing is manufactured of a structural steel having 
the following characteristic constants: 


pg = 0.285 Ib./in.?* é = O44 B:bn./P.f. 

mh = O62 Xx 167° ad = 00155 m.*/sec: 
B.t.u./in. sec. °F. 

a =75X10%°R-! E = 29 X 10° psi 

Consequently, 1 + 2(h./hy) = 5 (100) 


and Eq. (Zo 


becomes 
(101) 


f 


This means that the temperature of the juncture ¢ 
cover plate and web is 20 per cent lower than would be 
the temperature of the same cover plate under the 
same conditions of heating if the webs were absent. 
This result is correct only if the inequality, Eq. (2S), 1s 
satisfied. Substitution yields 


t < 87/16 XK 0.0155 = 258 sec. (102) 


if the criterion is applied to the web. The cover plate 
would yield a value 9/4 times greater. 
From Eq. (37) the average temperature of the cover 


plate is 


Tree = Ht — 0.00310") °F. (103) 

From Eq. (44) B=6 (104) 
and from Eq. (46) 

M* = 6/7 = 0.858 (105) 


Consequently, the average stress in the cover plate is 
on the basis of Eq. (47): 


Care = 30.8SHt(1 — 0.0219¢'*) psi (106) 


The uth harmonic of the thermal stress distribution 
is obtained from Eq. (62) as 


jn = —(1,910/n?) Ht ?@(A) (107) 


In the determination of the buckling conditions, the 
value of the reference stress oy defined in Eq. (68) is 
needed. It has the value 


oa) = 2,850 psi (108) 


TICAL SCIL 


NCES NOVEMBER, 1956 
Under uniformly distributed stresses, the panel of th, 
cover plate bounded by two shear webs will buck; 


when C = 4 as was shown in connection with Eq. (§ 
11,400 psi 109 


oO —= lo = 


The analysis presented in the body of the paper has 
shown that small deviations from uniformity in th 


stress distribution have a small effect on the averay 
compressive stress at buckling. Eq. (106) contains 
the information regarding the average stress. Th 


critical value, of course, may be reached after a short 
time of heating when // is large or after a long time oj 
heating when // is small. If it is known that the heat 
ing rate and the heat capacity of the cover plate ar 
such that the rise in temperature in the cover plate ir 
the absence of webs would be 


H = 40°F./sec. 110 


then, after the lapse of 10 sec., the temperature rise in 
the middle of the cover plate is about 400°F., and, from 
Eq. (106), 


On c= 11,480 psi 11] 


As at that time, Eq. (58) yields 


A = 0.203 112 

and, from Fig. 4, &(A) = 0.027 113 
one obtains from Eq. (107 

Jz = —1,682 psi ll4 


From Eq. (5), the eifect of the second harmonic on 
the buckling stress is expressed by the divisor 1 — 
(p/2) which in this case is 1.071. Hence, the actual 
buckling stress is 


o-, = 11,480/1.071 = 10,700 psi 115 


buckling would occur 
therefore after 9.33 sec. of heating. At that time, of 
course, the value of j2 is slightly different from the value 


In a second approximation, 


just calculated, but the determination of a third ap- 
proximate value is not warranted because the differ 
ence is very small. 


Substitution in Eq. (2) yields 


gq = 0.7170 B.t.u./ft.? sec. 116 
The value of // given earlier as 40°F. /sec. corresponds 
therefore to 

g = 28.68 B.t.u./ft.? sec. 117 


Under a constant heat input of this magnitude, the 
cover plate buckles after less than 10 sec. As this time 
is much shorter than the 25S sec. calculated in Eq 
(102), the simplifying assumptions made in the analysis 
regarding the length of the plate were justified. 
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Conical Techniques for Incompressible 
Nonviscous Flow 


HOMER J. STEWART? anp ALLEN I. ORMSBEE: 
California Institute of Technology 


SUMMARY 


41 analytical investigation was made of incompressible poten 
fields which are homogeneous of order zero. Expressions 
such fields were derived, rules for superposition of such 
re de veloped, and certain superpositions were made, re 
sulting in expressions for the flow fields associated with constant 
h source and vortex sheets of finite extent 
mstant strength source sheets were applied to the con 
f aerodynamic models of thin nonlifting wings of polyg 


The constant strength vortex 


uction ¢ 
lan form and cross section 
were used to construct models of thin wings with constant 


ressure distribution 


SYMBOLS 


real constants 


pressure coefficient 


2u/q 

= free-stream velocity 

= Var se +s 

complex velocity functions 

the real parts of l’, V, IV; perturbation velocities 
parallel to the x, y, and s axes, respectively 

= the imaginary parts of U’, V, IW, respectively 

= Cartesian coordinates in physical space 


= { Vy + 7 y + r 
= wi(x+r 
x+? 


v vector operator |(0/ON), (O0/O¥), (0/02 


INTRODUCTION 


._ STUDY is concerned with the properties of 
those solutions to the incompressible three-dimen- 
sional potential flow equations which are homogeneous 
of order zero. 

Donkin! has shown that the general homogeneous 
solution of order zero to the Laplace equation in three 


limensions may be written in the form 


S = F(e) + Gia) (1 


where e= (y + 22)/(% + Vx? + yy? + 2? 
éis the complex conjugate of «, and F and G are arbi- 
trary functions. 

Jones? and Lampert* have studied subsonic potential 


flow fields in which the velocity components (wu, 7, w) 
are homogeneous of order zero, and by superposing 
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certain of these fields were able to obtain flows satis 
fying certain of the boundary conditions for a finite 
body. Some of Lampert’s basic solutions were 1m- 
proper due to the presence of sidewash and or down- 
wash discontinuities which did not cancel out in the sub 
sequent superposition for finite bodies 

Flow fields in which the velocity components are 
homogeneous of order zero—so-called ‘“‘conical’’ flows 
are characterized chiefly by the appearance of sheet dis 
turbances source sheets or vortex sheets—which 
cover infinite planar sectors in space, and line disturb 
which 


€.g., 


ances—e.g., line sources or vortices have a 


linear variation of strength and extend to infinity 
Their 
the absence of a characteristic length permits no insight 


direct physical significance is questionable since 
into the nature of the flow at “‘large’’ distances from 
such disturbances. 

It can be shown, however, that under certain cir 
cumstances, proper superposition of such flows, in 
which a characteristic length is introduced, reduces the 
extent of the disturbances and allows a discussion of 
the flow at large distances from such disturbances. In 
these cases, the appropriate boundary conditions can be 
applied. 

In the present treatment, complex conical “velocity 
functions” (U’, V, W) are defined, the real parts of which 
represent the physical velocity components (“, 7, w 
Relations governing U, V’, and II’ are derived, restric- 
tions are discussed on the types of singularities these 
functions may have, and techniques are developed for 
the superposition of conical fields to obtain nonconical, 
physically significant flows. Several specific conical 
fields are discussed, and superpositions of these fields 
are made to provide examples of three-dimensional 


fields of interest to the aerodynamicist 


THE CONICAL VARIABLE € 
Before proceeding to a discussion of specific conical 
solutions to the subsonic flow equations, it will be of 
interest to examine certain features of the conical vari- 


able 


Referring to Fig. 1, we have 


II 


tan - O< Bsa 
>) 


and 


arg e = tan—'(2/y) = ¢ 0O<¢ 
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Thus, considered as a mapping of the xyz space onto 
the « = & + 7m plane, half-rays emanating from the 
origin in the xyz space are mapped onto points in the e 
plane. Furthermore, for each half-ray there exists a 
unique pair (8, ¢) for which O < 8 < wand0 < ¢< 27; 
and, consequently, since tan (8/2) is single-valued on 
this range, each half-ray is mapped onto a unique point 
in the « plane. The inverse transformation 


B = 2tan— |e <= 2 = 
g = arge CS e<s 


is one-valued by the restricted ranges on 8 and ¢ and, 
consequently, for every point in the ¢ plane (including 
the point at infinity) there is a corresponding half-ray 
in the xyz space. Therefore the transformation is one- 
to-one in both directions, in the sense that for each 
half-ray in xyz space there corresponds one and only 
one point in the ¢ plane, and conversely, for each point 
in the ¢ plane there corresponds one and only one half- 
ray in the xyz space. In particular, it is noted (see 


Fig. 2) that (a) The plane z = 0 is mapped onto the 
real axis of ¢, (b) The plane y = 0 is mapped onto the 
imaginary axis of e, and (c) The plane x = 0 is mapped 


onto the unit circle in the « plane. The one-to-one 
nature of the transformation allows, in the following, 
the treatment of the variable ¢ indiscriminately as 
representing either a point in the complex plane or a 
half-ray in the wyz space. 


THE COMPATIBILITY RELATIONS 


Suppose that U(e), V(e), W(e) are analytic functions 
of e, and that they are separated into their real and 
imaginary parts as 


U=ut in*, V=90 + w"*, W = w + iw* 
Then, by Eq. (1), 
V24= V0= Vw=0 
where 
V? = (07/dx?) + (07/dy?) + (07/027) 


It is asked under what conditions (wu, v, w) may be 
interpreted as the velocity components of a_three- 
dimensional incompressible flow field. A necessary and 
sufficient condition for this is that the complex vector 
F = (U, V, W) satisfy 


vx FF = (3) 


The real parts of these relations are the physical condi- 
tion of irrotationality and continuity. The necessity of 
Eqs. (2) and (3) is due to the assumption of analyticity 
on U, V, and W, and may be shown by a direct applica- 
tion of the Cauchy-Riemann equations in the ¢ plane 
to each of the complex functions L’, V, and W sepa- 
rately. 

Tt Eqs. (4), (5), and (6) may be derived as a by-product of this 


process. The details of this derivation are omitted since they are 
straightforward but algebraically quite complicated. 
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Since V?U = V?l’ = V?W = 0, by Eq. (1), and since 
these functions are homogeneous of order zero, Eq. (2 
is satisfied if Eq. (3) is. Considering Eq. (3) by com. 
ponents then, suppose that 


OU /dy = OV/ox 


Then 
(dU /de)(O€/Oy) = (dV /de)(O€/Ox) 
But 
O€/Ox = —(e/r) 
and 


Oc/Oy = (1 — e?)/2r 
Therefore, 
adU/de = [2 (6? — 1) |(dV/de) 4) 
Similarly, the conditions 


OU /dz = OW/ox Ol’/os = OW’ /dy 


imply 
dU/de = [2ie/(e2 + 1) ](dW/de) (5 
and 
dV/de = i[(e? — 1)/(e? + 1) ](dW de (6 
Since e(v, vy, 2) and its derivatives are single-valued, 


the above arguments may be inverted and Eqs. (4), (5), 
and (6) are a necessary and sufficient condition for the 
fulfillment of Eq. (2). It is to be noted that any one 
of Eqs. (4), (5), and (6) may be derived from the other 
two, so that this system is not determinant. These 
relations are referred to as the compatibility relations 
and are directly analogous to compatibility relations of 


supersonic conical flow. 


SOURCE FIELDS 


We now proceed to a discussion of a basic solution 
to Eq. (3) that may be interpreted as a particular dis- 
tribution of sources in the xy plane. This basic solu- 
tion is 


sss 
’ 2W Eo éE— € 
U,(e) = oe log 

w(é€o” + 1) e + (1/6) 


: woffeo” — | e—& | 
J 1(€) = . log — log € 
ri\ e+] e+ (1/6) f 


: —1Wy jle + (1 €o) (e — €)) 
Wile) = log - 
sia T . | € f 


where ¢€) and wy are positive real constants, ¢) <1. The 
associated real velocity field is 
2WoE ety 


& = Gi (0) = log 
we? +1) ~ le + (1/6) 


% = (VY; = 


Wo €0" — | ¢€¢=— £6 (7a 
( x ) log + log |e | 
Tv €0~ + l € + ( l €0) J 


Wo € + (1/€) J(e — €) 
Ww, = GR! (Wi) = arg 2 eo) I(e 


1 é 
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show 
spact 
value 
syste 
cuts 
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CONICAL 


In order to make w, single-valued, the ¢ plane is cut as 
shown in Fig. 3(a). The corresponding cuts in xyz 
space are shown in Fig. 4. 
valued, it is not necessary to construct a specific cut 


Since “, and 7 are single- 


system for LU’; and ;, however, a convenient system of 
cuts for these two functions is shown in Fig. 3(b). 
These cuts are not significant in the real velocity field 
and are deleted from Fig. 4 for clarity. 

If eis made to approach a point on the cut to the right 
of the origin from above, it is seen from the expression 


for w; in Eq. (7a) that 


lim w, = w if O<i< « 


e—i+10 
whereas, if approached from below, 


lim WwW) = — W if 0 «.. é . € 


e—t—i0 


The same values of w, occur on the cut to the left of 
—|/«@. Elsewhere on the & axis w; has the value zero. 

The field of Eq. (7a) is thus the field associated with a 
uniform source sheet of strength wy, distributed over a 
This field is represented 


double-sector in the xy plane.T 
schematically in Fig. 5, where the values of w are shown 
in the xy plane, approached from above. 

This particular field will be referred to as F(x, y, 2; 
«) Where F, is the vector with complex components 
(Ui, Vi, Wi). 

Formally, fields of the type Fi(x, y, 2; €) may be 
superimposed to produce fields of finite extent. Such 
a superposition is one which gives a triangular source- 


sheet and may be expressed as 


G(x, y, 3; €, €) = (1/2) [Filx, 9, 2; €0) — 
Fi(x —c, y — b, 2; «) + F(x —c, y — 6,251) — 
Fi(x — c, y, 2; 1)] (8) 


where c is a real positive constant and 


This superposition is shown graphically in Fig. 6, 
where the values of w in the xy plane, approached from 
above, are given for the component fields and the re- 
sultant field. The values of w on the xy plane ap- 
proached from below are the negative of those shown. 

It should be noted that 56 as defined above may be 
either positive or negative, depending on the magnitude 
of &, and in particular, —G,[x, y, 2; (1/€), ¢] is the re- 
flection of Gi(x, y, 2; €, ©) across the x axis. 

The question is raised at this point as to whether such 
a superposition is proper, since the function Gi(x, y, 2; 
, €) is not strictly defined along the half-rays asso- 
ciated with the singularities of the conical field F,. 
These half-rays are indicated by the dashed lines of 
Fig. 6(e).] Certain general statements will be made 
with respect to this question in the following section. 


t Jones, in reference 2, apparently did not consider explicitly 
the singularities in this field at — 1/¢ and infinity and in his initial 
superpositions assumed w = 0 far upstream—i.e.,¢ = ©. This 
difference cancels in superposing to obtain finite span and chord 


bodies and his solutions for these cases are correct 
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SUPERPOSITION OF CONICAL FIELDS 


The general requirement in superposing conical fields 
is the reduction of disturbances of infinite extent to 
disturbances of finite extent. The procedure, as indi 
cated in the formal superposition above, is to super- 
impose identical conical fields in pairs, one of the pair 
being displaced, in xyz space, from the other along a 
half-ray that is a singularity of the two fields. Thus 
the first two conical fields (a) and (b) of Fig. 6 are identi 
cal, field (b) being displaced from field (a) along the sing 
ular half-ray « = «, and the superposition is a subtrac 
tive one. 

Consider now a conical field F(x, y, 2) = (€), say, 
and suppose it is desired to superimpose this field on it 
self by the procedure indicated above. For convent 
ence, we shall suppose that the xyz axes are oriented 
so that the singular half-ray in which we are interested 
lies along the negative x axis, corresponding to 
The superposition may then be expressed as 


G(x, y, 2; c¢) = F(x, y, 2) — F(x — oc, y, 2 Q) 
or equivalently 
Pie, § lle — f(< 
where 
¢ = (Vy Tr 72 [fx —-c+ V(r -—cf?+ y+ 2 


and the notation g(e; ¢) is used to emphasize the fact 
that « and ¢ are not independent variables, their values 
being uniquely determined by the specification of the 
variables (x, y, z) and the parameter c. 

The geometry associated with this superposition is 
shown in Fig. 7, and it is seen that the half-rays e and 
¢ intersect at the point (x, y, 2), the « ray emanating 
from the origin (0, 0, 0) and the ¢ ray emanating from 
the point (c, 0, 0). Now suppose that the point (x, y, 
z) is not on the x axis and that the value of ¢ is fixed. 
If the point (x, y, ) is allowed to move out on the € ray, 
¢ will vary continuously, and it is clear from the geome- 
try that 

lim ¢ = « 


€ const 
If the point (x, y, s) is on the x axis, there are three possi- 
bilities: 


For y = 2 = 0 


fi) sx we have « = ¢ = 0 
iz; ec ser we have « = Oand¢ 
(3) x <0 we have e = ¢ = « 


It should be noted in connection with the third possibil- 
ity, however, that 
lim (e” — ¢") 
2-0 
0 
0 
does not exist. 
It is also of interest to note, in connection with the 


three possibilities above, that 
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CONICAL 


la mrs ct 
lim (€/¢) = (x% — c)/x 
z—0 
9a) forO<x <e 
lim (€/f) = 0 
a0 
ja forx <O 
lim (€/f) = x/(x — c) 


20 


Returning to the fields F and G of Eq. (S) it is seen, 
as a consequence of the geometrical properties discussed 


above, that if «is a point of regularity of f(e), then 


lim G(x, y,2;¢0) = lim [f(e) — f(g)] = 0 


— ti—e 


For the particular case where (x, y, 3) is on the x axis, 
the value of G is defined as 


G(x, 0,0;c) = lim (F(x, 9,2) — Fw —c, 9,2 


where this limit exists. 
It was assumed to begin with that F(x, y, 2) was 


singular on the negative x axis (e = ©); and since 
lim (e” — ¢ 
20 
x 0 
0 


does not exist, G(x, 0, 0; c) will not exist if the singular- 
ity of F on the negative x axis is in the nature of a pole 
or an isolated essential singularity of f(e) at « = 
If poles are present in the field, a more complex super- 
position is necessary. 

On the other hand, if f(€) has a branch point at e = 


such that f(©) can be made to exist by introducing 
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suitable cuts, or if f(€) has a logarithmic singularity at 
« = ©, then G(x, 0, 0; c) will exist for x < 0, provided 
that the proper cuts are introduced. 

The above discussion was carried out for the case of 
a singularity along the negative x axis; however, the 
results are quite general and may be applied to any 
other half-ray since any half-ray may be made to coin- 
cide with the negative x axis by a simple rotation of the 
space. 

It is to be noted that, in the example given above, G 
is singular on the segment 0 < x < ¢ which corresponds 
to « = 0 and ¢ = This result is general and the 
effect of superpositions of the type discussed is to cancel 
all but a finite segment of the singular half-ray of the 
conical field. 

Since in general a given conical field will have more 
than one singularity, several superpositions must be 
made in order to provide proper cancellation of all sin- 
gularities; however, the linear nature of the process 
permits the superposition to be carried out in any order. 

As an example, consider the field F\(x, y, 2; €)) of Eq 
(7) and the superposition indicated by Eq. (S) and Fig 
6. Since the field PF, possesses only logarithmic singular 
ities, it qualifies for superposition by the procedure out 
lined above. Referring to Fig. 6, it is seen that the 
combination of (a) and (b) provides proper cancella- 
tion along the singularity « = e, (a) and (d) provides 
proper cancellation along the x axis, (b) and (¢) pro- 


7 


vides proper cancellation along the line y = 3, and (c) 
and (d) provides proper cancellation along the line x = 
c. Thus, the resultant field (e) is nonsingular except 
on the boundaries of the triangle, and the resultant 
velocities (u, v, w) become arbitrarily small at large 


distances from the triangle. 


NONLIFTING THIN WINGS 


The conical field /i(xv, y, 2; €), Eq. (7), and super 
positions of the type Gi(x, y, 2; €, ¢) of Eq. (8) may be 
used to construct fields associated with a particular 
class of physically interesting bodies. This class con- 
sists of symmetrical wings at zero angle of attack whose 
maximum thickness is small enough that the boundary 
conditions on the surface of the wing may be replaced 
by boundary conditions in the chord-plane of the wing. 

Although an infinite superposition, or integration, of 
such fields may be used, in principle, to obtain the field 
associated with a thin wing of arbitrary plan form and 
thickness distribution, the technique is particularly 
suitable for determining the flow about thin wings of 
polygonal plan form and airfoil section, such as are 
common in much of present-day supersonic wing de- 
sign. The subsonic characteristics of such wings are 
not irrelevant, since invariably these wings must accel- 
erate through the subsonic regime before attaining 
supersonic speeds. 

As an example, consider the field 


G3(X, Y, 3° €0) C1, Co) = (1/2) [Fi(x, y, 2; «€ — 
J { 3 


Fi(x —c a b, 2; € — Fi(x = Co. 7, 3, 6 — 
1 ly 2 J 
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CONICAL 


where b [2eo/(1 — e07) Ja 


ind ¢, €2, and €) are real positive constants. 

This field, illustrated in Fig. 8, represents the flow 
past a swept thin wedge of finite chord and span, fol- 
lowed by a constant thickness infinite chord plate. 

The reflection of this field across the x axis is given by 
(1/€0), C1, C2] 


—G3[x, ¥, 3; 


and the superp¢ sition 


(a(X, V, 2; € ( ( ‘ G3(X, Vy 2; €0, C1, C2 ss 

GsiX, V, 2; | € » C1, Co] — Gs(X — Ce, ¥, 8; €0, Ci, Ce) 
. , oF - - | 
G3[x — C3, y, 2; (1 /€0), C1, Ce] 


provides the perturbation field for a swept, constant 
chord wing with a modified double-wedge airfoil sec- 
tion (see Fig. 9). The pressure at any point on the 
surface may be calculated by a superposition of values 
of u, computed from Eq. (7a). 

Values for the pressure coefficient along the mid- 
chord line and along the half-panel chord line are plot- 
tedin Fig. 10 for the wing above, for the case c 34/3, 
V3. This corresponds to 
1 sweep angle of 30° and an aspect ratio of 6. 


Co 3+/3/2)c3, and e l 


VORTEX FIELDS-—CONSTANT PRESSURE WINGS 


We now consider conical fields which may be inter- 
preted as lifting surfaces. 

As before, the boundary conditions are linearized. 
That is, the angle of attack is presumed to be small and 
the pressure perturbations are linear in the streamwise 
perturbation velocity u.f Because of these approxi- 
mations, the perturbation velocity “ must be an odd 
function of z, as shown by general linearized potential 
wing theory, and « must be zero on s = 0 except on the 
wing itself, where it will be of the same magnitude on 
the top of the wing and on the bottom, but of opposite 
sign. A further consequence of linearization is that the 
vorticity vector must lie in the xy plane and, since u 


0) off the wing on z 0, 


Ov/Ox 0 on z = 0 

The sidewash v in the wake is independent of x. If, 
however, v is also conical, then v must be constant in 
the wake. 

As a consequence, the field associated with an infinite 
lifting sector possessing both a leading and trailing 
edge is not conical, since the induced velocities as de- 
fined by the law of Biot and Savart do not exist for a 
constant-strength infinite vortex sheet. At- 
tempts to describe such a field as conical lead to singu- 


sector 


larities elsewhere in the field. 

T In this discussion the wing will always lie “in” the xy plane 
ind the streamwise direction will be coincident with the positive 
Y AXIS 

{ Lampert, in reference 3, satisfies the boundary conditions on 
(for such a field, however singularities in v associated with his 
solution violate the requirement of zero vorticity except on the 
wing and wake, and these singularities do not cancel in his super- 


positions for finite wings 


T 
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Consider now the conical field F2(x, y, 2; €9) defined 


by 
Us (lUo/ 7) log [(e eg eT & 
V2 — (1/27 [(e 7-— | log ee T 2 log «| 
Ws — (uy /27e0) [(eo” + 1) log (e® — €”) — 2 log €] 


1, % 1s a real posi 
to the 


where ¢€ is a real constant, 0 < « 
tive constant, and the e plane is cut from « € 
left through real negative values of «€ to « 

This field was discussed by Lampert and the real 
part of this field consists of vortex sheets in the xy plane 
as shown in Fig. 11, where the values of u and v shown 
are those on the xy plane approached from above. 


The superposition 
G;(x, ¥, 2; €0, C) F2(x, y, 2; €0) — Fe(x — c, y, 2; €) 


where c is real and positive yields the perturbation field 
for an infinite span swept wing with constant pressure 
everywhere on the surface (see Fig. 12). The reader is 
referred to Lampert’s paper for a discussion of the cam 
ber distribution associated with this wing. 

In order to construct finite-span wings, it is neces 


sary to consider the field F3(x, y, 2; €) defined by 


Us = [iuo(1 — e)/m]flog [(e + 1)/(e — «)] 4 
{(1 + €&)/(1l — « | log i(e€ ] e— ¢ ' 
Vs — [2u0(1 — €07)/mwe | log [e/(e — « 


7 Uo(1 — €0*) € € 
Ws tos ) t ( x 
Te €—e l — « 


e— ] € € 
log + log ( 
et l — <€ € 


0< ea <1 
where the ¢ plane is cut on the real axis between « ] 
and e« —1. 


This field consists of vortex sheets in the xy plane 
as shown in Fig. 13. The reflection of this field across 


the x axis is given by 
— F3(x, y, 2; —« 

The superposition 
Ge(x, ¥, 2; €0, C) F;(x, y, 2; €0) — 
where b = 200 


and c is a real positive constant, gives the field indi- 
cated in Fig. 14. The reflection of this field across the 


¥ axis is given by 
—Ge(x, ¥, 2; —€0, ¢ 
Now consider the superposition defined by 


2)Ge(x, ¥, 2; €, C) — 


G7(x, y, 2; €, c) = (1 ‘ 
€0G5(X, ¥, 2; 


(1/2)Ge(x, v, 2; —€0, C) 


LL @ 


This field is shown in Fig. 15 and is the perturbation 
field for a constant pressure delta wing. The surface 


shape necessary to support this pressure distribution is 
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a highly cambered one and is characterized by logarith- 
mic infinities in slope at the leading and trailing edges 
and along the centerline. 

The reverse-delta constant pressure wing is obtained 
by the superposition 


Gs(x, ¥, 2; €0, C) —(1/2)Ge(x, y + 5, 2; €, ¢) + 


(1/2)Ge(x, y — b, 2; —e€0, c) — e0Gs(x, y, 2; 1, €) 


where 0 is given in Eq. (11). This field is indicated in 
Fig. 16. 

Since any polygonal plan form which does not have 
streamwise tips (parallel to the x» axis) and which is 
symmetric about a streamwise axis can be partitioned 
into isosceles triangles like those of Figs. 15 and 16, 
these fields may be combined to obtain a constant pres 
sure wing of such a plan form. In such a superposi- 
tion, the singularities in downwash in the component 
flows will cancel except on the boundaries of and down- 


stream of the vertices of the resultant wing. 


CONCLUDING REMARKS 


The foregoing represents a very convenient method 
of treating polygonal source sheets and constant pres- 
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The 
transformation may be used to extend this treatment 


sure regions in potential flow. Prandtl-Glauert 
to subsonic compressible fields. 

In addition it would appear that the constant pres 
sure delta wings might be superposed to provide an 
approximate model for the general lifting wing. Pre 
liminary attempts to construct such a model, using a 
scheme similar to the horseshoe vortex superposition 
of Falkner,* were not successful. The difficulties en- 
countered appeared to be chiefly due to the presence 
of singularities in downwash associated with the edges 
of the triangular vorticity patches and the lack of a 
rational basis for locating downwash control points 
Work in this direction is being continued. 


REFERENCES 


1 Donkin, W. F., Phil. Trans. 1857. 
Bateman, H., Partial Differential Eg 
Physics, Dover Publications, New York, 1944 

2 Jones, R. T., Subsonic Flow over Thin Oblique Airfoils at Zer 
Lift, NACA TN 1340, 1947 

’ Lampert, Seymour, Conical Flox 


Referred to on p. 357 of 


uations of Mathematu 


Methods Applied to Un 

formly Loaded Wings in Subsonic Flow, Journal of the Aeronaut! 
18, No. 2, p. 107, February, 1951 

M., The Solution of Lifting Plane Problems by 

Theory, Aero. Research Council, R&M No. 2591 


cal Sciences, Vol 
‘ Falkner, V. 
Vortex Lattice 


1953 


On Turbulent Flow Near a Wall 


(Continued from page 1011) 


REFERENCES 


1 Stokes, G. G., On the Effect of the Internal Friction of Fluids 
on the Motion of Pendulums, Trans. Cambridge Philos. Soc., 
Vol. 9, 1851. 

2 Laufer, John, The Structure of Turbulence in Fully Developed 
Pipe Flow, NACA TN 2954, June, 1953. 

3 Nikuradse, J., Gesetamdssigkeiten der Turbulenten Strémung 
in Glatten Rohren, VDI, Forschungsheft 356, 1932 


‘ Nikuradse, J., rauhen Rohren, VDI, 
Forschungsheft 361, 1933 

5 Taylor, G. I., Conditions at the Surface of a Hot Body Expose 
to the Wind, Tech. Rep., Advisory Committee for Aeronautics 
Vol. 2, R&M 272, May, 1916 

6 Prandtl, L., Bemerkung iiber den 
Physik. Zeitschr., Vol. 29, 1928 

7yvon Karman, Th., The Analogy Between Fluid Friction ane 
Heat Transfer, Trans. ASME, Vol. 61, 1939 


Strémungsgesetze in 


Warmetibergang im Rohr 





sree eo 


plan 
ones § 
of lift 
Incre; 
ducin 
incre¢ 
Rec 
1496 


* Ac 


Glauert 
-atment 


it pres- 
vide an 

Pre- 
using a 
position 
Hes en- 
resence 
e edges 
‘k of a 
points 


or Af 
. dos Ol 


ematica 
at Zer 


to Um 


ronaut! 


VDI, 


‘x pose d 


autics 
Rohr, 


m and 








Reduction of Drag Due to Lift in Supersonic 
Flight by Distributing Lift Along a Fuselage 


R. M. LICHER* 
Douglas Aircraft Company 


SUMMARY 
It is shown that it is possible to reduce the drag due to lift for 
n isolated wing in supersonic flight by adding a fuselage carrying 
» proper distribution of lift; this is possible even if the wing 
yy itself has the optimum distribution of lift for its plan form 
The optimum fuselage lift distribution for any given wing lift 
jistribution is given along with simple upper and lower bounds 
n the possible drag reductions. One method of calculating the 
lrag of the wing-fusclage combination is outlined 
An elliptic wing carrying the optimum lift distribution for that 
plan form is used to demonstrate the drag reduction possible 
These reductions can be quite large—of the order of 30 to 40 
per cent in some cases. It must be emphasized that it may not 
be possible to generate the required fuselage lift distribution with- 
ut producing excessive vortex drag, and that experiments are 
necessary to determine how much of the theoretical improvements 


can actually be achieved in practice 
(1) SYMBOLS 


wing span 
= wing chord 


Dware = Wave drag due to lift 
D(6) = wave drag of an equivalent lineal lift distribution 
dy = fuselage length 
(0) = length of equivalent lineal lift distribution 
E = complete elliptic integral of the second kind 
= subscript referring to fuselage 
| = integral defined in Eq. (23 
K = complete elliptic integral of the first kind 


= total lift 


fx 
bal --|: 
0 


In (1 — ¢)dt, Euler’s dilogarithm 


equivalent lineal lift distribution 
average equivalent lineal lift distribution 
VW = Mach Number 


= dynamic pressure 


subscript referring to wing 
= coordinate in free-stream direction 
£ = Jacobian zeta function 
ps =«7V/6-1 
9 = polar coordinate measuring angular position from 
horizontal on Hayes control surface 
mn = Bb/¢ 


(II) INTRODUCTION 


| gears OF STUDIES have dealt with the problem of 
reducing drag due to lift on various isolated wing 
plan forms. For certain wings, such as the elliptical 
ones studied by R. T. Jones,! the optimum distribution 
of lift to produce minimum drag due to lift is known. 
Increasing the wing span is known to be effective in re- 
ducing the vortex drag; in a somewhat similar manner, 
increasing the wing chord can help to reduce wave drag 
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due tolift. But because of other considerations, such as 
skin friction and structural limitations, it is not always 
possible to use as large a wing as might be desired for 
reduction of drag due to lift. Hence it is of interest to 
see how the effective chord of the whole airplane can be 
increased by distributing lift along the fuselage. 

One might consider the plan form of the wing-fuselage 
configuration to be given and ask what is the optimum 
distribution of lift on both wing and fuselage for mini- 
mum drag. A somewhat easier but more restrictive 
problem is to fix, in addition to plan form, the wing lift 
distribution and ask for the optimum fuselage lift dis- 
tribution. The latter problem is studied here. 

Linearized supersonic flow theory is used throughout 
this report and the drag calculations are based on the 
method of W. Hayes.” 
interference drag formulas developed by E. W. 
For simplicity, the fuselage is represented by a 


Extensive use is made of the 
Gra- 
ham. 
line in the plane of the wing parallel to the free stream. 
The net lift carried on the ‘‘fuselage”’ line must be zero; 
otherwise the vortex drag would be infinite. For the 
drag calculations in this study, it is assumed that it is 
possible to load the fuselage in the optimum way. 
From the practical viewpoint, however, it is not known 
how much lift a fuselage can be designed to carry ef- 
ficiently. Undoubtedly some lift can be distributed 
along the fuselage, but perhaps the full amount calcu- 
lated in this report cannot be realized without exces- 
sive vortex shedding and consequent increase of vortex 
drag. However, these calculated drag reductions are 
large enough in some cases so that the attainment of a 
fraction of this reduction might be important. 

FOR A 


(III) Optimum FUSELAGE LIFT DISTRIBUTION 


GIVEN WING LIFT DISTRIBUTION 


Once the wing lift distribution is given, the optimum 
fuselage lift distribution can be found relatively easily. 
It is somewhat more difficult, however, to calculate the 
drag of the optimum combination. The development 
of the criterion for the optimum fuselage lift distribu- 
tion given below is a specialization of the proof given by 
M.E. Graham.‘ She considered the general optimiza- 
tion problem when both lift and thickness distributions 
in space are involved. For the planar case considered 
in this report there is no interference between lift and 
thickness distributions, so only the lift case will be de- 
veloped here. 

The evaluation of wave drag is based on the method 
of Hayes’ in which the drag is related to the momentum 
outflow through a cylindrical control surface at a large 
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distance from and surrounding the given configuration. surface. d is the length over which /(x, 6) exists: jj the « 
The lifting configuration, viewed from each angle 6 the lineal loadings are considered to lie along the fuse. In 
on the cylindrical control surface, can be represented lage axis, d may be either the fuselage length or the term 
by an equivalent lineal lift distribution (lift distributed length of the wing Mach envelope depending on the wave 
along a line parallel to the free stream). Hayes has Mach number and geometry (see Fig. 1). Hayes has wing 
shown that the drag of the actual configuration is equal shown that the /(x, #) are determined by passing through net Ii 
to a weighted average of the drags due to lift of the the configuration a series of parallel planes inclined at to lil 
equivalent distributions. Thus, the Mach angle with respect to the x axis and lumping duce: 
1 all lift intensities intercepted by any one plane at the 
(2 | D(6) sin? 6 dé point where that plane cuts an arbitrary axis parallel IV 
— a (1) to the X-axis. Te 
Die) = — B’ | | ) (E, a)I"(x, 8) X rhe lineal distribution Ux, 6) can be separated inti On 
l6mgd —a/2 9 the part ly, due to the wing alone, and /,, the fuselage rn 
ta |x — Eldtdx lift distribution which is independent of @: “ 7 
where /’(x, 0) = di/(x, 0) /dx and /(x, 6) is the equivalent I(x, 0) = lw(x, 6) + I(x) D) fopt\" 
lineal lift distribution as a function of the streamwise the | 
coordinate x and the angular position @ on the control In Eq. (1) the @ integration can be carried out first. — 
lusela 
9 > | ever 
2n %lT 2r 
[ I'(x, O)1'(E, 0) sin? 6. dé = | l’w(x, 0)l'w(&, 0) sin? 6 dé + 1',(E)l’ (x) [ sin? 6 dé + : gat 
70 e /7 0 | ind 
I’ (&) / l'w(x, 0) sin? 6.dé + 1';(x) , l’w(&, 0) sin? 6 dé | Optin 
‘ ‘ 
Define the following average wing lineal lift distribution: | a 
gn 
F l [- a s intens 
lw(x) = lw(x, 6) sin? 6 dé 3 
TJ O that 1 
Then plan f 
) 9 due te 
Qn Qn 
i} I(x, 0)1'(E, 8) sin? 6. dé = [ l'w(x, O)l'w(&, 6) sin? 6 dO + ml’ El’ (x) + ml’ (EN w(x) + al’ p(x) w(8) tusela, 
0 /0 In thi 
2x the wi 
= { l'w(x, 0)l'w(&, 0) sin? 6 d0 — al’w(x)l’w(é) + wf[l'(&) + ptin 
0 I 
I'w(€) [0 (xe) + Uw must I 
Eq. (1) can then be written 
_g (27 a2 paz — 
Dae | f tl’ w(x, Ol’ wl, 0) — l'w(x)l'w(E) + 
16mr*q/ 0 —d/2J —d/n 
[l’w(x) + lox) [2 w(E) + (8) ]f In wx — —) dxedé sin? @ dé (4 
; 
ee a Fn 
The first term in Eq. (4) is the wave drag of the wing the average wing distribution. Since the first tw 
alone; the second term is the wave drag of a lineal lift terms are specified for the optimization problem con- 
distribution, which is just the average wing lift distri- sidered in this report, the wave drag will be a mini- | Where 
bution; and the third term is the wave drag of a lineal lift mum when the third term is a minimum. Thus the 
distribution, which is composed of the fuselage plus optimum fuselage lift distribution is one which forms 
in combination with the average wing distribution, an The no 
optimum lineal lift distribution. If the fuselage length Fried 
is greater than the length of the wing Mach envelope The of 
then the average wing plus fuselage distribution must is the ¢ 
be an optimum for the given fuselage length, d,, and 
total wing lift, Lw. It is known that the optimum 
lineal lift distribution under these conditions is elliptic. Thus tl 
However, if the average wing equivalent distribution 





extends beyond the fuselage length, then the optimum 


average wing plus fuselage distribution must meet the This di 
added constraint of having its shape prescribed over This fu 





part of its length; this is a more difficult problem, but drag is 











REDUCTION OF DRAG 1039 
ists; if the optimum distribution generally can be found. 
‘ ; . a4 , wu z 
1e fuse. In general, it should be possible to make the third ae 4 jini anaest Bidiaccn 
or the term of Eq. (4) smaller than the second so that the | 
on t -e drag of the configuration is less than that of the 
ne wave 5 | 
ves has wing alone. If the fuselage is required to carry zero 
hrough net lift, so that its vortex drag is zero, then the drag due 
ined at to lift (the sum of wave and vortex drag) can be re- 
um ping duced by carrying lift on the fuselage. 
at the 
arallel ’) EXAMPLE ILLUSTRATING POSSIBLE REDUCTION Fic. 2 
| 
IN DRAG DUE TO LIF1 
ed into ) the lift distribution is given over a wing plan 
Mace the ) n over ;¢ s pl: = lB oh tt teh AE A as 
uselage wate See oe SI Phe equivalent lineal lift distribution for the elliptic 
form, its equivalent lineal distributions, /j(x, 4), can ': : a . 
Les} ; ' ; ' wing when viewed from any angle 6 on the Haves evlin 
be calculated; so also can the average, /w(x), an . . . = ; 
oo a abies 5 Ai drical control surface is elliptic: 
, x), the optimum fuselage lift distribution. Once 
5 he loadings are known, there are several methods inh 3 9 
ead loadi gs oO ire eve ul met 10ds lw(x. @) = 16Ly red A ]* DP sae [2x d(6 5) 
= est which may be used to calculate the drag of the wing- 
irst. : — é 
selace ruration. ‘ ‘ acticz “ases, } 4 ) 
fuselage ¢ nfigurati n ; For most practical cases, how — d(0) = cV/1 +(@b/c)? cos? 8 6) 
ever, the calculations will probably not be easy to carry 
out analytically. One example that can be completed d(@) is the length of the wing equivalent lineal loading at 
ynalytically is discussed in this section. any angle 6. The average wing equivalent lineal lift 
= ; ; sh se ad * distribution is, from Eq. (3), 
in” 6 dé Optimum Fuselage Loading for Elliptic Plan-Form Wing 
Consider an elliptic plan-form wing in supersonic l6Ly € 1 + w* cos? 6 — 4& sin? 6dé ei 
: . . lw(x) = (7) 
fight of span 6 and chord c, which carries a constant W\*) ate Je 1 + uw? cos? 0 ; 
J intensity of lift (see Fig. 2). R. T. Jones! has shown 
that this lift distribution is optimum for the isolated where § = x‘ 
plan form. It is of interest to find how much the drag w= Bbc 
cae ‘ ‘ . - >) a 
due to lift can be reduced by distributing lift along a . ju- ee 
/ "= , aa a 
£) iuselage of length d; > cV 1 + (8b/c)?, 8 =V M2 — 1. leos™!'V (4? — 1) /u 2&3 
> ° ° ° e 
In this case, the fuselage is longer than the length of ( Vi1l+ 2p? 
the wing average equivalent lift distribution; thus, the 
optimum fuselage plus average wing lift distribution lw(x) extends over the length d(@),,,, = ¢V 1 +(8b c)? 
I’ w(x must be elliptic. Eq. (7) can be evaluated in terms of elliptic integrals as 
l6LwV 1+ pw? — 42 V a? — y*K(a)Z(6, a) l 
- K(a) — — Ea) <= ¢ =< . 
a ss * a — 2 
7 Cu YV1— ¥’ 2 
lo (4 l6Lw Vil+p?— 42) . /1 aV a® — y? K(1/a@)Z[n,(1/a) | we Sane 
: af’ k — —ak ) <§@<s-V14+,? 
ron a yV1l-y a ys 2 
(8) 
st two 
m con- 
. mini- Where y = p/V1 4+ p? a=pV1itpw — 42 drag, the interference drag between wing and fuselage 
us the : ss Matt is negative and the total drag of the configuration will 
: 6 = 1 1 (4 ) =s 1 (a) 7 c 
Pe — ue 7 —_ } be shown to be less than that of the wing alone. 
ion, afl lhe notation for the elliptic integrals is that of Byrd and “ - 
Rie ] . . “er ay Drag of timum Configuration 
length Friedman.’ /y(x) is plotted in Fig. 3(b) for » = 1. ow fig 
velope The optimum distribution of lift Lw over a length d, The wave drag of the optimum lifting elliptic wing 
1 must is the elliptic one and fuselage has been shown to be equal to the wave 
l,, and drag of three separate configurations: 
timum lop (Xx) = (4Lw/mds)V 1 — (2x/d,)* (9) 7 a 
— T 4 “ sae . : , - Doan = Doing wave Tt D, ng+fus. wate ~—— D, ng wave (4) 
Lliptic. hus the optimum fuselage lift distribution is 
bution The first term, the wing alone wave drag, is known:! 
semen I(x) = lope) — lw(x) (10) 
Thie dicts: . : , i , Duing wave = (Lw?/aqb*)(V 1 + pw? — 1 (11) 
et the This distribution is shown in Fig. 3(c) for d,/c = 3. sh an gb*) | . 
1 over This fuselage loading carries zero net lift so its vortex The second term is also known, since the fuselage lift 
m, but drag is zero. Although the fuselage has some wave distribution has been chosen to make the average wing 
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C) OPTIMUM FUSELAGE L/FT OISTRIBUTION 
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Fic. 3 


plus fuselage distribution an optimum for the given 
fuselage length and net lift carried. For this example, 
the optimum distribution is elliptic and its drag is 

= B*Ly? 
Diving + fus = 


(12) 


ar 2rqd/ 

The only remaining drag term to be calculated is 
that of a lineal lift distribution which is just the average 
wing lift distribution lw(x), Eq. (8). The 
method of superposition of optimum shapes developed 
by E. W. Graham® is particularly suited to the drag 
evaluation of this example because of the form of /j(v). 
It has been pointed out that /;(x, 6), as a function of x, 
is an elliptic distribution for each angle @—that is, it is 
an optimum lineal distribution of lift. /j(«) is just a 
weighted average with respect to @ of /w(x, 0), [Eq. (3) ] 
so can be considered as being composed of an infinite 
nested set of optimum elliptic distributions of lift. By 


—-1.€., 


B? J Ly? 


hina wate = 
oe 2mq \d(0)? 


Integrating by parts and noting that 


28°Lw? 
Deine wate — : = 
1Wwqc~ JO 
9 9 
D = 9P 48°Lw* 
wing wate ~— —4wing wave ae 
1W°qc~ 
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nested it is meant that each elliptic distribution ejthe, 
completely contains or is completely contained by a; 
other. The drag of each individual lift distributioy ; 
known and, for such a nested set, the interference dra 
expressions have a simple form. 

For a nested set of 7 lineal elliptic lift distributions 
each of which carries lift L; over a length d,, d; < 
the drag is shown in reference 3 to be 

n } 1 ? 
D = >> L/Dop(d;) +2 dD L; Dd LiDopld 
1 ] j+1 
where D,,,(d;) is the drag of a unit net lift distribute 
in the optimum manner (elliptically) over the lengt! 
d;. The explicit form of the drag expression for th 
elliptic distribution is 
Dopd ) = 8? 


2mrqd;* 


Subtistuting this in Eq. (13) and rewriting the summa 





tions in a form more suitable for this example (as is 
done in reference 6), 


D = (82/2xq){[(9.L,)2/d,?] + 


? 1 j 
> (> L,)? (1/42) — (/djy1”)]] 15) | 
] l ; 


In reference 6 it is shown how Eq. (15) may be extende 
to include an infinite number of nested optimums by 
taking the limit as 1; becomes infinitesimally small an 
djs41 > d; + Ad, The length of each distribution 1 
the set is [Eq. (6) ] 


d(g) = cV 1+ nw’ cos’? ¢ 


is a parameter which varies from 0 to 7? 
of the distributions is 


(4Ly 


where ¢ 
The lift carried in each 


Tio) = 7) sin? ¢ 


The total lift carried in the nested set is 





L = in yon eee 
Noting that 
DL = [" Lea | 
and that 
(1/d;?) — (1/dj41?) > d(1/d?) = —2[d(g) ]~*(dd/de)d¢ 


then Eq. (15) can be written as 


"Title f*" *- —2(dd/dy)d¢ 
= [ | J] sin? xd || —_ lt (16) | 
JV T Vo da? 
r/2 sin? g dg 
1 + up? cos? ¢ 
(17 


r/2 (9 — — 
f (2¢ — sin 2¢) sin’ gd¢ 
0 1 + yp’ cos? ¢ 
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Phe integral in Eq. (17) can be evaluated as 


TI, = {[20. + w?)Le)/u2} — { [C1 + uw?) In ( 


where 
- Adé6 l r: a a ae | =a 
Ji + y2cos?6) 2V1+ 42 2 pe AW + Vi + ntl I 


—f1 In (1 — ¢) dé 


Lol Xx = 
For x < 1, the following series representation of {o(«) may be used :’ 
+. oe 
L(x) = Dn *x ex 4 


If the resultant expression for the wave drag of the average wing is substituted in Eq. (4), the wave drag of the con- 


figuration can be expressed as 


ipa = 


~~ 


— D + (46°L yy ™q / 20 


The elliptic wing vortex drag is 
D . re = Ly 1qb* (?] 


As a measure of the reduction in total drag due to lift (wave plus vortex) possible through distributing lift along 
the fuselage, the ratio of wave plus vortex drag of the configuration to that of the wing alone is given Eq. 22 and 
plotted versus the reduced aspect ratio, 6:R, in Fig. 4 for various values of the fuselage length to wing chord ratio 


/ 





(Dy P Dw) wav eon SAS YD = 7) 2a d,)? —~VI1 + pw? + (4/2°) — 
[4(1 + pw?) In (1 + w2)/e2u?] + [SC + w2)Jo/m?]} (22) 
For 6.R = O the wing vortex drag dominates (wing Slender body theory can be used to get a preliminary 


wave drag goes to zero) so that the curves all start 
from unity. In Fig. 5 the results are replotted versus 
Mach Number for three wing aspect ratios. 

Eq. (22) is valid only if the fuselage is longer than the 
wing Mach envelope (d;/c > ‘1 +2); thus some of the 
curves of Figs. 4 and 5 cover only a limited Mach Num- 
ber range. The results could be extended to higher 
Mach Numbers but the computations would be more 
difficult because of the more complicated form of the 
optimum average-wing-plus-fuselage lift distribution. 
For V1 + pu? > d,/c, 
the fuselage becomes less able to effectively modify the 


as the Mach Number increases 


average wing lift distribution and so the curves of Figs. 
tand 5 will eventually approach unity again. Thus, 
lor a given configuration, it appears that the fuselage 
can be used most effectively in reducing drag due to 
lift when the Mach Number is such that the length of 
the wing Mach envelope is somewhat less than the fuse- 
lage length. 

In Eq. (4) the wave drag expression is broken into 
three terms, two of which are fixed when the wing lift 
distribution is fixed; the third term is the drag of the 
This last 
distribution may always be made optimum if the aver- 


average wing plus fuselage lift distribution. 


age wing distribution is contained within the fuselage 
length. Thus the elliptic wing may be shifted fore and 
aft on the fuselage to some extent without changing 
the drag values given in Fig. 4 (although the required 
luselage lift distribution will change). 


estimate of the camber necessary to produce the re- 
quired lift on the fuselage. As an example, choose the 
following parameters: J/ = 2, :R = 1.47, andd,/c = 3. 
It can be seen from Fig. 4 that the drag due to lift could 
be reduced 32 per cent if the optimum fuselage lift dis- 
The required lift dis- 
The total positive lift 
carried on the front part of the fuselage can be found 
If the 


fuselage is assumed to have a circular cross section and 


tribution could be developed. 
tribution is given in Fig. 6. 


by graphical integration to be AL = 0.22 Ly. 


a fineness ratio of 15, slender body theory requires 
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that the air stream must be deflected by an angle a = 
ho Sera Thus 


a C, of 0.2 would require a change in local angle of at- 
The length of fuselage 


degrees to produce the required AL. 


tack over the fuselage of 36° 
carrying the positive load is about 97 per cent of the 
wing chord. The required angle is, of course, much too 
high to be practical and is not within the range of angles 
in which slender body theory is valid. Hence some 
device for making the fuselage a better lifting surface is 
needed. If, for a given volume, the fuselage cross sec- 
tion were made elliptical rather than circular, an in- 
crease in lifting surface could be obtained. Doubling 
the fuselage width would cut the required local angle 
of attack down by a factor of four. 

Another possibility is to accept a smaller drag saving. 
If a nonoptimum fuselage lift distribution is used in 
order to cut the positive AZ in half (AL = 0.11 Ly), 
the resultant drag increase over the optimum case will 
be less than 10 per cent of the wing drag; thus a saving 
in drag due to lift of more than 22 per cent still results. 
Cutting AZ in half also halves the required change in 
angle of attack, according to slender body theory. If 
now the fuselage cross section is changed to a 2:1 ellipse 
while maintaining the same fuselage cross-sectional 
areas, the required local angle of attack is reduced to 
9°; this value is probably within the range of accept- 
able values for practical configurations. 

It is still necessary to consider the effect of vortex 
shedding from the fuselage. In the calculation of drag 
it has been assumed that no net lift is carried on the 
fuselage, and that no vorticity from the fuselage ap- 
pears in the Trefftz plane. Actually vortices may be 
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shed from different parts of the fuselage and those repre- 
senting negative lift may be displaced vertically from 
those representing an equal positive lift. This would 
produce a vortex pattern in the Trefftz plane which 
has no vertical momentum (no net lift) but which has 
energy, corresponding to drag. Such a drag penalty for 
carrying a lift distribution on the fuselage could seri- 
ously decrease the overall drag savings obtainable. 
Experimental data are necessary to determine how 
efficiently the fuselage can carry lift, especially at 
angles of attack where vortex shedding may be impor- 
tant. 
boundary-layer control to prevent vortex separation. 


It might be necessary to consider the use of 
Once the actual fuselage effectiveness is determined 
experiments with a properly designed wing-fuselage 
would be desirable to determine the actual drag reduc- 
tion that can be obtained. 

(V) Some SIMPLE UPPER AND LOWER BOUNDS FOR 
DRAG DUE TO LIFT FOR OTHER CONFIGURATIONS 


As mentioned in Section (III), it is usually easier to 
determine the optimum fuselage lift distribution for a 
given wing than it is to calculate the drag of the optt- 
mum configuration. Hence it is desirable to have some 
simple lower and upper bounds for the optimum drag 
values that can be used to estimate the possible im- 
provement. 

Probably the simplest lower bound is the drag value 
which would be achieved by assuming the wave drag to 
be that of the fuselage alone if it carried all the wing 
lift in an optimum manner [distributed elliptically, see 
Eq. (12)]. Adding to this the wing’s vortex drag gives 
the following result for the lower bound to the drag of 
the lifting wing-fuselage combination: 
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REDUCTION 


dD, B= dD, ng vorter + (@°Lw? 2arqd/ (23) 


For the upper bound, consider a wing of arbitrary 
plan form carrying lift Lw and a fuselage of length d, 
Let d,, be the length of the Mach 
Each wing equivalent 


as shown in Fig. 7. 
envelope enclosing the wing. 
lineal lift distribution will be enclosed within this 
length d,,, which will also be the length of the average 
wing equivalent lineal distribution. Pick for the fuse- 
lage lift distribution the difference of two elliptic dis- 
tributions of length d, and d,,, each of which contains a 
total lift L,; thus the net fuselage lift is zero. The 
interference drag formulas of reference 3 can be ap- 
plied to this nested set; they hold even if the innermost 
lineal distribution (the wing in this case) is not an op- 
timum. Defining Dw(Lw, @) as the wave drag of each 
wing equivalent lineal lift distribution and D,(L,, x) as 
the wave drag of an optimum distribution of lift L, 
over the length x, the upper bound for the wave drag at 
each angle @ can be written as [see Eq. (15) ] 


DALw, d;) + DA(Lw — L,, d,,) — 
DALy L,, dr) + Dw(Lw, 6) — D,(Lw, d») 


Di®up 


From Eq. (1), the wave drag is 


Disa weve = } D(@) sin? 6d6 
0 


= DALy , d,) + DLi “a Ly, d, ws 
DALy — = d;) T Proins ware 
DA Lw, d») (24) 


The form of D; is given in Eq. (12); adding in the wing 


vortex drag 


dD, a = Li *D | l, d,) + (Lw — L,)?[D,(1, dm) 
D,(1, d;)] + Dwing. — Lw?D,(1, dm) (25) 


If Eq. (25) is differentiated with respect to L, and the 
result set equal to zero, the value of L, for the mini- 
mum upper bound of this form can be found; it is 
easily shown that 

Ly Of as Ly 
so that 


D, a, = D, = DLw, d») + DLy, dy) (26) 
= Dying — (B?Lw?/22q)[(1/dm?) — (1/d,’)] 


These simple bounds, Eqs. (23) and (26), are plotted 
in Fig. 8 versus reduced aspect ratio for the elliptic wing 
along with the curve of Fig. 4 for d,/c = 4. An im- 
proved but more complicated upper bound taken from 
reference 4 is also plotted in Fig. 8. The reader is refer- 
red to reference 4 for a more thorough discussion of 
upper and lower bounds for this type of problem. 


(VI) CONCLUSIONS 


It has been shown that in theory the drag due to lift 
lor a given wing in supersonic flight generally can be 
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Fic. 8. Upper and lower bounds for reduction of drag due to lift 
on a fuselage plus elliptic plan-form wing 


reduced by carrying a distribution of lift along a fuse- 
lage. This is demonstrated for elliptic plan-form 
wings where as much as a 30 or 40 per cent reduction in 
drag is found for some cases. 

However, it is suggested that the optimum lift dis- 
tributions required cannot, in practical cases, be carried 
efficiently on a fuselage. By cambering the fuselage 
and by making use of the change in the fuselage cross- 
sectional area, and perhaps through other means, it 
should be possible to develop some lift on the fuselage. 
However, the possibility of vortex shedding from the 
lifting fuselage, with its attendant increase in vortex 
drag, must also be considered, unless some means, such 
as boundary-layer control, is developed to prevent it. 
Even if only a fraction of the improvement in drag 
noted can be achieved in a practical case it might prove 
important. It would be of interest to perform experi- 
ments to determine the amount of lift intensity a fuse- 
lage can develop efficiently and how nearly the optimum 
fuselage lift distribution can be approached. 
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The Aerodynamic Force on an Airfoil in a 
Moving Gust’ 


JOHN W. MILES7 


Unwersity of California, Los Angeles 


SUMMARY 


The increment of lift due to a sharp edged gust that moves 
across a two-dimensional airfoil at a relative velocity different 
from the flight velocity is calculated on the assumption of incom- 
the results to helicopter 


pressible flow. The application of 


blade-wing interference is discussed 


(1) INTRODUCTION 


- apsiese AN AIRFOIL that moves with a velocity 
(’ and is subjected to an upwash field moving 
with velocity (1 — A~')U so that the latter moves 
downstream or upstream as J is posi- 
'U. Fur- 


along the airfoil 
tive or negative with a relative velocity \ 
ther, let x be a dimensionless, streamwise coordinate 
having the values —1 and +1 at the leading and trail- 
ing edges, respectively, and o a dimensionless time 
variable defined according to 

o = (velocity:time) semichord = Ut/b (1.1) 


Then the upwash along the airfoil may be placed in the 


form 

v(x, o) = Uf (o — Aw — |A}) (1.2) 
where / is the instantaneous upwash angle—on the 
usual assumption |v) < l~—defined such that the 


upwash meets the leading (v = —1) or trailing (x = 
+1) edge at time zero, as X is positive or negative, 
respectively (see Figs. la, 1b). It is required to deter- 
mine the resulting lift on the airfoil as a function of the 
dimensionless time o and the relative velocity param- 
eter X. 

The problem posed in the preceding paragraph re- 
duces to the well-known, sharp edged gust problem if 
\ = | and to the indicial lift (‘Wagner’) problem if 
\ = 0, both special cases for which complete solutions 
are available.’ ? The more general problem, \ + 0 
or 1, may arise in treating blade-wing interference on a 
winged helicopter {see Section (S)] or in determining 
the forces on a turbine rotor blade moving through the 
wake of a stator blade. In the latter connection, Kemp? 
already has obtained the required solution when the 
prescribed upwash v In the former 
connection, the aperiodic or transient problem also is 
of interest and will be considered in the following dis- 
Insofar as only the total lift on the airfoil 


is sinusoidal. 


cussion. 
is required—rather than the pressure distribution 

Received November 4, 1955. 
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airfoil (A > 0). 


the result for an arbitrary time dependence might by 


calculated directly from Kemp's result by Fourier 
superposition; but, in order to exhibit more clearly the 
separation of response into ‘‘circulation”’ and ‘‘virtual 
mass'’ components, it is expedient to start from 
Schwarz’ basic result’ for the pressure on an oscillating 
airfoil having an arbitrary, chordwise distribution of 


upwash. 


(2) THe OSCILLATING AIRFOI 


Suppose that the prescribed upwash v has the har- 
monic time dependence exp (wf) where w is the angular 
frequency. Converting to the dimensionless time o, as 
defined by Eq. (1.1), the exponential appears as ex 
(tka) where k is the reduced frequency defined by 


U 21 


Then, in terms of v(v, o) and k, the pressure jump 


across the wing, as given by Schwarz,’ reads 


pl’ ("* : 
Ap = P| y=0- — P| y=04+ = | J>1C(k) (l + 
7 7 { 
7] OL(A, o) 
cos ¢) — cos ¢] tan + + 
2 Od 
ik sin d L(6, ot v(cos¢,a)do@ (2.2 
where the trigonometric change of variable 
x = cos 9 2.3 
the ‘‘Theodorsen function” 


C(k) = [Ah (k) + tH (k)J-! INO (Rk) (24 
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Sharp edged gust moving downstream relative 
Fic. 1b. 
relative to airfoil (A < 0). 
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Sharp edged gust moving upstream 
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and the logarithmic influence function 


TI — cos (0+ ¢)- 
1(@. o) = In | = 
” 1 — cos (0 — ¢) 
a: > 


— 
, ] 


n~! sin (76) sin (n@ 


have been introduced. 


A modified form of Schwarz’ result that proves more 


convenient in many applications follows from Eq. (2.2 
after integrating by parts the term v(OL O¢); the result 


may be placed in the form 


2pl (5) is 
Ap = tan - | 


"le ee 
‘(cos @, a) dé + [ sin @ L(6, o) X 
Tr J0 


[C(k) (1 + cos ¢) — cos ¢| X 


Dv (cos ¢, o) F 
do 
Dt 


(2.6 


where Dv Dt represents the ‘‘total acceleration,’ viz. 


Di _ Ov Ov U fov . U Di 
l +r 7 + tkv)}) = (2. 
Dt O(bx) Of b \oOx b De 


[he first term in Eq. (2.6) has a distribution over the 
airfoil that is identical with the pressure distribution 
over a flat plate in steady flow tan (6 2) = 


1 — x) (1 + x)] “—exhibiting an integrable singu- 


note: 


larity at the leading edge, x = —1, and satisfying the 


Kutta condition of vanishing pressure perturbation at 
the trailing edge, x = +1. The second term in Eq. 
2.6), depending on the total acceleration distribution, 
might be designated as the contribution of rate of change 
of virtual momentum to the pressure, but, as brought 
out below, this definition is to some extent arbitrary. 
Independent of this definition, it should be remarked 
that (Dv Do) vanishes for a stationary upwash pat- 
tern—i.e., one that depends only on o — x [A = 1 in 
Section (1)] rather than on o and x separately. 

The lift coefficient corresponding to Eq. (2.2) or 


Eq. (2.6) is given by 


l 1 Ap 
(1/9 — d(bx) = 
l 


2b. p 2] 
T Ap ; 
— sin 6d6 (2.8 
J0 (p l a 
and takes the alternative forms 
2 fF 
7 U | | Cie) (1 + cos @) v(cos ¢, a) + 
J0 L 
sin? } : (cos ¢, | dod (2.9a) 
Og 


— o 
= U ) [C(e) (1 + cos @)—cos ¢]-v(cos o, 0) + 


Dv 


sin? @ D (cos ¢, os dg (2.9b) 
Co 


Where ikv has been replaced by (0v/0c), and (Dzv/Do) 
is defined by Eq. (2.7). 
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The last terms in Eqs. (2.9a) and (2.9b) will be 
designated, respectively, as the ‘‘partial’’ and ‘‘total”’ 
While this 
breakdown is largely a matter of convenience, it may 
be established that the lift coefficient calculated from 
on the assumption of a purely potential 


virtual momentum components of the lift. 


the upwash v 


flow—i.e., no circulation—is given by 


») 


Cc, = : / sin? @ (COS @, o do 
U Jo Oo 


(2.10) 


Accordingly, the breakdown of Eq. (2.9a) gains a cer- 


tain logical preference; nevertheless, as in the following 
results for \ > 0, it is not necessarily the more expedient. 


LAPLACE TRANSFORM SOLUTION 


(3) 


Let f(s) denote the Laplace transform of f(¢) with 
respect to the dimensionless time o defined according to 


[ e "fla)de (33. 
J 0 


1 f+, t; 
: deal f(s)ds ( 
271 « —% 


spectral 


f(s) = [f(e) = la) 


3.1b) 


flo) = L 1f(s) = 


In this transformation, the dimensionless 


parameter s appears as a generalization of the reduced 


twl 


frequency, such that s = 7k for real s—thus e“ = 
e’ = e*’. Accordingly, given the complex ampli- 
tude—i.e., the coefficient of e’“°—of the lift coefficient 
say ¢, (ik)—resulting from an upwash distribution of 
complex amplitude #(x, 7k), the Laplace transform of 
the lift coefficient due to an upwash distribution* 
having the Laplace transform (x, s) is given by ¢;, (s). 
Turning now to the specific upwash distribution of 
Eq. (1.2), the required transforms are given by 
Ue (i! +495 5) (3,2) 


OLX, S = 


we 
~ 


Oc = Si(x, S ee 


L Ov(x, o 


— A)si(x, S) (3.4) 


L Do(d, ¢)/Do = (1 


Substituting Eqs. (3.2)—(3.4) in Eq. (2.9) and replacing 
k by —is in the argument of the Theodorsen function 


then yields 


€1(s) 


= 2e~ '*'5 f(s) [C(—is) (1 + cos ¢) + 
70 


cos @ 


s sin? gle * do (3.5% 


= 2e7 X'S F(s) } }[C(—is) (1 + cos¢) — 
7 0 
cos¢?]+ (1 — A)s sin? gfe “8% dd = (3.5b) 


In order to reduce further the results of Eq. (3.5), it 
is expedient to rewrite them as follows: 


* Note that with respect to an observer on the airfoil the in- 
itial value of this upwash must be reckoned as zero. This as- 
sumption is implicit in the Laplace transform operation of Eq 
However, this assumption does not prevent the treat- 


(3.3). 
ment of such limiting cases as the sharp edged gust 
the footnote following Eq. (3.12 


See also 
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E(s) = 2m [ho(s, X) + hi(s, d)] sf(s) 


= 2m } [ho(s, A) + AAi(s, A)] + 
(1 — hi (s, r)} sf(s) 


(3.6a) 


(3.6b) 
= 27s h(s, d) f(s) (3.6c) 


where /) and /, are given by 


r RR ; l . a 
ho(s, X) = se '*'§ C(—i4s)- { grant 
TJ0 


(1 + cos d)dd@ (3.7) 
Als l ths ASCOS@ .°_9 — 
. e sin? @d¢d (3.8) 

/7 0 


us 


hy(s, A) =e 


The integrals in Eqs. (3.7) and (3.8) may be evaluated 
in terms of modified Bessel functions of the first kind 
{reference 5, 3.71(9)]; moreover, C(—7is), as defined 
by Eq. (2.4), may be re-expressed in terms of modified 
Bessel functions of the second kind according to 


(3.9a) 


C(—1ts) = [Ko(s) + K,i(s)]—! Ki(s) 


= 1 — Gs) (3.9b) 
where @(s) is the Heaviside transform—i.e., the La- 


place transform multiplied by s—of the ‘‘Wagner lift- 


deficiency function’ in the notation of Sears.” Thus, 
Eqs. (3.7) and (3.8) go over to 
ho(s, 3) = s—e7'*"§ [Iy(As) — 

T,(As)] [Ko(s) + Ai(s)]~! Ay(s) (3.10) 
hi(s, X) = (As)~!e7 "5 Tas) (3.11) 


The quantity / in Eq. (3.6) is the Laplace transform 
of the lift growth function in the important special case 
where f(c) is a unit step function—so that f = s~!. 
In this special case, the final lift coefficient is 27, and 
the lift growth function has been normalized by re- 
moval of this factor. Moreover, on the basis of the 
foregoing discussion, /y(o, \)—the inverse transform 
of ho(s, \)—may be termed the 
function” and f/,(o, \) the ‘virtual mass function.” 


The remainder of this paper will be devoted to the 


“circulation growth 


determination of fo and /, since the more general result 
may be recovered by Duhamel superposition; thus, the 
inverse transform of Eq. (3.6c) is given by* 


Ci(o) 


= = f(0 + )h(o, X¥) + 


eo 
| h(o — a1, A) f’(o;)do, (3.12) 
/ O04 


(4) THE VIRTUAL Mass FUNCTION 


The inverse transform of /,, as given by Eq. (3.11), 
can be evaluated exactly (reference 6, p. 124) to obtain 


* The result of Eq. (3.12) provides for a finite initial value of 
f(a) whereas this initial value has been assumed to vanish in 
Eqs. (3.2)-(3.4); thus, Eq. (3.12) might better be regarded as a 
limiting case by allowing f(a) to exhibit a steep rise from O to 
f(O0 +) in a very small time interval starting from a = 0. 
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Fic. 2. wiA\h(o, dX) vs. |A|~'o, as given by Eq. (4.2 


hy(o, X¥) = (wd?)— (21A « ~— 6)" O<a<2IA\ (41a 


lA 


= ( se 


The absolute value of \ must be used in this result jp 
consequence of the fact that 4; is an even function oj 
A. A more symmetrical result for /; is given by 


w!A|Ay(o, A) = 


o 
ioe 
r 


which plots as a semicircle, as shown in Fig. 2 


It is of particular interest to note that the total time 


integral of /, is independent of X, viz., 


eZ iA | 
hi(o, \) do = - 13 
JV 


Thus, /; consists of a semielliptical pulse of width 
!, and: total area 1/2; 


2)X|, maximum height (7A) ) 
in particular, it reduces to a delta function in the limit 
A= 8 viz... 
; # 
lim /,(o, \) = , 5(a) 1.4 
A\— 0 — 
This last term usually is omitted in plotting the lit 
coefficient due to a suddenly imposed unit angle of 
attack (the ‘‘Wagner problem’’), but this omission is to 
some extent misleading. 
POWERS OI 


(5) EXPANSION OF Jt) IN ASCENDING 


The “circulation function’? iyo may be adequatel) 
approximated by an expansion in ascending powers 
of « for small o and by a series of exponentials for 
large o. It is found that interpolation between these 
two approximations is required in the range A < ¢< 
2|\| —a check being furnished by the known results! ° 
1; for positive \ the combination hy + Mr 

gives a much 


for \ = 
following the breakdown of Eq. (3.6b) 
smoother fairing than does jy alone, whereas the con- 
verse holds for negative X. 

In order to expand J in ascending powers of ¢, al 
expansion of jp in descending powers of s is required. 
Introducing the asymptotic expansions of the modified 
Bessel functions [reference 5, 7.23(1)—(3)], it is found 
that 
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} 3 } SA S 2d 2567 
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— — + gs? + Q(A s—*) 5.la 
64 O4A 1024," S*r* 
0 \ ( \ 1/2 lit i(1 l ) ' ( l 3 
\ : / nay = Ti / ) S — § _ " > a 
2) s\ * 4fal ~ 162 
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+ _ + os Sa se D. 
64 64) 512? 512)A . — 
7 1/2 3/2 l I 3 
1>0: hols, X) + AMi(s, A) = (27d) “9 ~~ Li = s~lh + a S 
q SA 16 128A" 
( l 3 15 ) | : 
_ A. af. § + O(A7~* 574 2.2 
32Xr 12SA° 512d J 


inverse transforms of Eqs. (5.1la), (5.1b), and 


( + nd Jo 
2r 


Taking the 


1>0: Apolo, A) = 


\>0: Iolo, X) + Ali(o, A) = 


It should be emphasized that these results, having 
been derived from asymptotic series, are not necessarily 
convergent and, therefore, are not reliable beyond the 


point at which successive terms begin to grow. How- 


ever, they appear to be adequate in the range 0 < o < 
\ ; in particular, for \ = 1 the sum of the results of 
Eq. (5.3a) and of Eq. (4.1) is found to agree with 
Sears’ more exact result? to better than 1 per cent in 


a a“ 


eS 


(6) EXPONENTIAL APPROXIMATION FOR /ip 


If ¢ is not small, a useful approximation to the in- 
verse transform of s~! C(—is) has been given by Jones’ 


in the exponential form 


[-! )S 1 [Ko(s as Ki(s) ] i K,(s){ = 
1 — 0.165e 00450 0.335e 0.3000 (6.1) 
which corresponds to 
|K ‘s) K(s) ] 1 Ki(s) = si — 
0.165(s + 0.045) — 0.335(s + 0.300) —! (6.2) 


Moreover, since J,,(As) is analytic in the entire s plane, 
the Heaviside expansion theorem yields 

1) \ 
L )e T,(As) (s + @) lt = 


[,(—dXa) e *? u(lo — \A\) (6.3) 


Where u is the unit step function. Thus, assuming 


thato > \ , making use of Eq. (6.2), and noting that 


(2) | | 
l— o 
T 12d 


(5.2) then yields 
l 3 
——{1- .” 
70 tr 
l 5 
] —_ — 
560 iA 


ie) ) 
oe: 
320! 
15 1515 
112A $48 


)o +- O(A~? o* 
)o +- ((A~* o* 


5.3b) 


a 


a l— - o- — 
Gon \ =) / 


l 3 
( - , 
120d 1 


I) and J, are odd and even functions, respectively, the 


inverse transform of Eq. (3.10) is found to be 


ho(o, X) = 1 — 0.165 [Jo(0.045A) + 
0.045(¢—|4!) _ 9.335 [Zo(0.300A) + 


a> id 
1,(0.045A) ] e 


1,(0.300A) } e 0.30014 A (6.4 

For \ = 1 the approximation of Eq. (6.4) is found to 
agree with Sears’ more exact calculation’ to better 
than 1 per cent for o 2 (in which range /; = 0); 


accordingly, it seems reasonable to expect a similar 
when \ # 1. Since the results 
of the preceding section are not expected to be accurate 
for o appreciably greater than ||, the gap |A 

o < 2\\| must be filled by (graphical) interpolation. 
In order to bring out the nature of the interpolation 
[cf., first paragraph in Section (5) ], Ao, i, Ao + Ma (for 
\ > 0 only), and ty) + My are plotted separately as 
These curves show 


accuracy for ¢ > 2A 


functions of o/ | | in Figs. 3a—3h. 


hy as given by Eq. (5.3) inO0 < ¢/|\! < 1 and by Eq. 
(6.4) in o/|X 2, h, as given by Eq. (4.1) in 0 < 
o/|| < 2, and ty + Ay as given by Eq. (5.4) in 
0 < ¢«/A| <1, all drawn as solid lines; while the 
interpolations of fy and /iy + Ak, in 1 < a@/\A} < 
In addition, io + /; is plotted 


versus o in Figs. 4a, 4b (for X = 0, 4) = Oin o > 0). 


2 are 


drawn as dashed lines. 


In interpreting the results of Figs. 3 and 4, represent- 
ing the ratio of the instantaneous lift coefficient to its 
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final value of 27 during the entry of a sharp edged 
gust, it should be remarked that the high peak in / in 
the neighborhood of o = |X| is associated with the 
discontinuity in velocity at the gust front. 
practical applications, however, the gust shape, as 
in particular, as 


In most 


represented by f(a), will be smooth 
assumed in Section 3, the initial value will be zero; 
accordingly, the corresponding lift coefficient, as deter- 
mined by the Duhamel superposition formula [Eq. 
3.12)], may be expected to exhibit a peak that is either 
much less pronounced or entirely absent. 

In connection with the actual evaluation of the 
Duhamel integral, it should be emphasized that hy + 
Ah, and (1 — XA), should be treated separately for 
\ > 0; similarly, 4o and /, should be treated separately 
ford < 0. or may 
be so approximated 
the Duhamel integral involving / exactly, but the Du- 
hamel integral involving hy + Ah; (A > O) or ho (A < 0) 
generally will have to be carried out numerically or 
graphically. Alternatively, it might be feasible to 
take, at least approximately, the Laplace transform of 
/, substitute in Eq. (3.6), and invert the overall products 
by the methods applied to jy and /; herein. 


If f(c) is given in analytical form 
it may be possible to evaluate 


7) EXPANSION IN DESCENDING POWERS OF (¢ — A) 


While the approximation of the preceding section 
should be adequate for large values of a, it is of some 
interest to examine the expansion in descending powers 


ol ¢ or, more conveniently, ¢ — | A Introducing the 
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expansions of the modified Bessel functions in ascending 


powers of s [reference 5, 3.7(2), 3.71(15)], it is found 
that Eq. 


. r "a 
Ao(s, X¥) =e » 3 * E + In (5) af. ; in | is 
s\? ] s l 
1 + In = a y + In eo rh s + 


Of a In ’ Be A‘s- 


3.10) is approximated by 


where y is Euler’s constant. The inverse transform of 


Eq. (7.1) is given by 


l1—(o— \A))7' + 
} (2 — (1/2)A] — 21n [2(¢ — |v )]EX 
ao — |A!)-? +0 [(o — [A|)—* In? (o — JA!)] (7.2) 


hoo, %) = 


Since, by hypothesis, ¢ is large compared with A, /y 
vanishes identically in the approximation under con- 
sideration, and h = hy. 

It is evident from Eq. (7.2) that / is relatively insensi 
tive to changes of \ for ¢ > |X|, and this is borne out 
by Fig. 4. Physically, this result may be ascribed 
to the fact that the gust front is well removed from the 
airfoil whence the lift on the airfoil is almost independ- 
ent of the relative velocity of this front. 


(S) APPLICATION TO WINGED HELICOPTER 


A problem to which the foregoing results are appli- 
cable is presented by the repeated loading of a helicopter 
blade due to the upwash field of a wing directly under 
the blade. If the blade chord is sufficiently small 
compared with the wing chord, the wing upwash field 
may be calculated as if the blade were absent. More- 
over, if the blade span is large compared with the wing 
chord, the flow at any blade section may be assumed 
to be two-dimensional. Then, as illustrated in Fig. 5, 
if I’ is the forward flight velocity of the wing and 2 
the angular velocity of the blade, the rearward* flight 
velocity of a retreating blade element at a distance 7 
from the shaft when the blade is transverse to the 
flight direction—i.e., directly over the wing, assuming 


the latter to be unswept—is given by 


U = Qr — V = V[(r/r.) — 1] (8.1) 


and the parameter X, as defined in Section (1), becomes 


A=1— (V/Qr) = 1 — (r,/r (8.2) 

where , is defined by 
? = V Q) (S.3) 
As r is greater or less than 7,—corresponding to rear- 
ward or forward flight of the blade element relative 
to still air—\ is positive or negative, respectively. 
In this respect, it should be remarked that the small 


* In the problem that motivated the present study, the blade 
shaft was located near the wing tip, and the blade elements were 
e.g., blade mounted on starboard wing 
In other 


retreating over the wing 
tip rotated counterclockwise when viewed from above 
configurations, r would have to be reckoned negative over the 


advancing blade. 
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Fic. 5. <A blade of chord 24 retreating with relative velocity 
Qr over a wing of chord 2/,. 
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perturbation theory, on which the result of Eq. (2.2) is 
based, breaks down in the neighborhood of r = +r, or 
A = O since the gust velocity in this neighborhood may 
not be assumed small compared with the flight velocity 
(U’) of the blade. It also may be noted that \ ap- 
proaches — © as r approaches zero whereas the numeri- 
cal results of the preceding sections extend only to \ = 
—1, although the analytical results are valid in —A = 
—1. Neither of these shortcomings is regarded as seri- 
ous in the calculation of blade bending moments; in 
the first instance, the integration with respect to 7 
would smooth out the error near r = r,—just as the 
leading-edge singularity in thin airfoil theory is 
smoothed out by integration; secondly, it probably 
would be sufficiently accurate to neglect the interval 
0<r< (1/2)r. (A = —latr = (1/2)r,) in the bending 
moment integral, especially since r, is usually small in 
those problems of principal interest. 

Turning to the calculation of the dimensionless 
upwash f(c), as introduced in Eq. (1.2), it may be 
supposed that the true upwash v due to the presence 
of the wing is prescribed relative to the wing velocity 
V as a function of the ratio of the true wing coordinate 
x, (see Fig. 5) to the wing semichord )j, viz.,* 


v= Vii (xX) by) (S. } ) 
Now, from Fig. 5 it is seen that 


x, = Ort — blx + (AL/A] + m1 


0 - 
(S.0a) 

* The referee has emphasized that present knowledge of blade 
loading on an isolated rotor often may not be adequate for the 


prediction of the distribution function fj. 
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where x;°° is the value of x; upstream of which v may 
be neglected; introducing ¢ from Eq. (1.1), 


) 


a = b[(Qr/U) ao — (AX /)] +m" — be (85h 


Introducing Eq. (8.5b) in Eq. (8.4) and equating the 
result to Eq. (1.2) then yields 


; V (; ) (F d 
Ko —- AX — |Al) = =h -o= 4 
l | b l r ) 


Expressing the velocity ratios in terms of \ from Eqs 
(S.1) and (8.2), Eq. (S.6) goes over to 


| [ b 
Ko — Ax — |A = (< - ALG) x 


or, equivalently, 


l Che xX} , 
Hla) = — | fy a SS 
i \ rd, b, 
which completes the required determination of f(¢). 
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On the Damping in Pitch of a Helicopter 
Rotor in Forward Flight 


David F. Gebhard and Leonard Goland 

Research Associate and Chief Project Engineer, Respectively 
Forrestal Research Center, Princeton University, Princeton, N.J 

June 7, 1956 


SYMBOLS 


blade section lift curve slope, 1 /rad 
1 blade coning angle, rad 
angle between perpendicular to tip path plane and control axis in 


longitudinal plane; positive rearward, rad 


angle between control axis and thrust vector; positive rearward 
rad 
b number of blades 
t angle between perpendicular to tip path plane and control axis in 


lateral plane; positive to right, rad 


B blade tip loss factor 


CH rotor force coefficient perpendicular to control axi positive 
rearward 

C7 thrust coefficient 

h distance from rotor hub to fuselage center of gravity, ft 

I blade flapping moment of inertia, slug ft 

VW pitching moment about the center of gravity; positive nose up, ft 
Ibs 

7 pitching velocity; positive nose up, rad. /sec 

R blade radius, ft 

S rotor disc area, ft 


rotor thrust; positive up, lbs 
blade mass constant 


average blade drag coefficient 


A average collective pitch angle, rad 

N inflow ratio, total flow parallel to control ax QR; positive when 
flow is down through disc 

“u advance ratio, flight speed/QR 

f air density, slugs. /ft 

o rotor solidity 

Q shaft angular velocity, rad. /sec 


bye HANDLING QUALITIES of present-day helicopters are sig 
nificantly governed by the damping in pitch of the vehicle 
a major role in producing this damping, and 
contribu 


The rotor plays 
consequently the accurate determination of the rotor 
tion is essential. This contribution is composed of the moments 


due to the resultant thrust inclination (or rotor drag) and flapping 


hinge offset. The thrust moment, to which this paper is devoted 


may be written as 


OM/d0qg = hT(0a‘’/a l 


In reference 1, the above factor, 0a’/Og, has been presented to 
equal 


) 


0a’/Og = — (27/y2){1.0 0.29 (@a0/C 2 


It was stated therein that this expression is fairly accurate for 
both pitch and roll for values of advance ratio up to w = 0.5 
However, using the above expression, difficulties were encoun 
tered by the authors in correlating flight test data for a Sikorsky 
S-55 helicopter with theoretical analysis, which indicated that a 
more detailed evaluation of the effect of advance ratio on this 
rotor damping term was necessary The results of this investiga 


tion are as follows 
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It is therefore suggested that Eq. (13) be used for those cay — 
er in which rotor thrust damping plays a significant role ae 
” . 
> REFERENCES (3) € 
 -4000 3 Or 
x Amer, Kenneth B., Theory of Helicopter Damping in Pitcl R : 
“ | EQ.I3 Comparison with Flight Measurements, NACA TN 2136 October, 1950 iy 
8 = 2 Gessow, Alfred, and Myers, Garry C., Jr., Aerodynamics cf the Heli F 
i lst Ed., The Macmillan Company, New York, 1952 . 
| ed 
ve | 
» —2000 —+ | 
be 
a 
Ww ae 
= -—!000 t t T | , . : a/f 
re) Application of Second-Order Shock-Expansion 
a Theory to Several Types of Bodies of 
4 O Revolution ( 
= 0 os 410 AS 20 DB DBD 3s 
FE ADVANCE RATIO ~ uw Robert E. Lavender and Raymond A. Deep 
Technical Feasibility Studies Office, R & D Division, OML se 
FIGURE |. VARIATION OF PITCH DAMPING WITH Redstone Arsenal, Huntsville, Ala ; 
ADVANCE RATIO June 18, 1956 
The rotor damping moment about the fuselage center of SUMMARY 
gravity may be expressed as Second-order shock-expansion theory is utilized to obtain equations for J 
é -e Cc » slope yitc y oO 2 irve slope { c 
OM og = hapoeS(QR)?] (0(2Cy/ac) da,](a1/dq) + the initial normal force curve slope, initial pitching moment curve sl 
‘or ; ; aor ) = and zero-lift wave drag for several type bodies of revolution. Bodies cor tribut 
0(2CH/ac)/Ob;|(0b Og) + [OA2CH/aa)/Oqg); (3 sidered are the cone-cylinder, cone-cylinder-frustum, cone-cylinder-frustum 
T , ~ . ; yooster, cone-frus ‘ -one-frus booster 1 
The rotor drag force coefficient, according to the method of DeSsine, GANS Themtees, ane cube aetem Sows 
reference “, page 196, and including a blade tip loss factor and a fisreObUCTION — 
term involving pitching rate, may be written as : wy allies 
: _ BODIES considered are shown in Fig. 1 and have a con On 
» we (ks 2/9 if 3/2) — NS ay oh. = . in — 
<Ca/ee (6uB — (6a, B* , (ud9B i : diameter of unity. Coefficients presented are based upor 
(3\a,B?/4) + (uai?B?/4) — (aob,B*/6) + the cross-sectional area of the cone and cone diameter as refer 
2 2 es ¢ e _— —s 
(a0%4B?/4) (qaoB*/62) (4 ence area and length, respectively The initial center of pressur 
Consequently the derivatives may be evaluated as in cone diameters from the nose is given by 
O(2Cy ag)/Oa; = (aB3 3) + (3\B? 4) + (ua, B? 4 (5 bs = | —dCy da (dCy da 
yi ) = — 3/R) . 
O0(2CH/aa)/ob; (ao B3/6 (6 DISCUSSION 
O(2CH/aa)/Oy = — (a)B3/6Q) (7 (1) Cone 
The blade flapping derivatives with respect to pitching velocity Values of pressure ratio, P;/P , Mach Number, ’ and initia 
are normal force curve slope, (dCy/da)., are obtained from the » oO 
Sane nal : | (5) © 
a 1¢ Massachusetts Institute of Technology cone tables rhe pite i On 
a ) ° P ‘. s - ‘ ' I 
=— (8 ing moment, including the effect of axial force, is given by 2 
Og BtyQ(1 — p?/2B? , WwW 
db 1 (dCy da). = —(2L./3) (dCy/da), sec? 0 
=— g -Cyli 
dg OL + p2/2B? ( (2) Cone sennanal 
On a cylinder following a cone: 
The steady-state coning angle, a), and the longitudinal flapping ‘ : Kil 
9 (dCy/da)., = (A) K,) (1 — e*'% - j 
angle, a;, may be expressed as (see reference 2, page 194 7 ‘ 
P - , (dCy da) = —(A,/K;) Pa + (1/KA, = 
ayo = (yB*/2)|(0/4)(1 + w?/B2) + (A\/3B) (10 y . s ' 
L. + L., + (1/Ki)Ie . 
u|(80/3) + (2d/B)] 
a= : (11 where A, = 4(A2/A1)At 1) | 
Bil — (u?/2B?)] 
. Set : : : (P/Po)/OS| » . es < 
Substituting the above expressions in Eq. (3), the expression for a Ts , P./P O1 
° -— (P, ) | 
the rotor thrust damping moment becomes dane | d 
aM hapaS(QR)? § 16 The functions necessary for evaluation of the constants A; an 
i apo S(QR)? 
== a 
oq 2 (25B4|1 — (u2/2B?2)] 
(CF 3\B2 am) paoB 12 ———_—_—————_— —1 vhere 
3 a 2 T W201 + 2/9 B2 ia i“ | 
3 4 2 120[{1 + (u?/2B?)] f AG, | 
or, using a tip loss value B = 0.97 (similar to reference 1) and aw Rey? 0.5 
i . : | a-<it Ge n ! ee" —- i = at | 
dropping small terms, >» XK 1,2 3,4 5,6 | 
OM/dq = — (1Qbh/R)(2.750 + 6.388 + k-~- L, —* ba me i ap = > . 
4.25ua,)(1 + 0.53u2) (13 All 
whicl 
The variation in pitch damping with advance ratio for the pre- - 
viously mentioned helicopter is presented in Fig. 1 based upon - : eas 
. ‘ a —_— L ' mM 
Eq. (13) and reference 1 The effects of reverse flow were not — t OF Re an 
included in deriving either expression, hence it was not deemed <6c | 
. ° 9 e <4 ——-+—- ——_ — 1 
advisable to carry the calculations beyond uw = 0.3. It will be >X 17 7" 
noted that there is a significant increment between the curve ob- a © . ‘ a! Lik 
. . ° . . 2 c f se 
tained using reference 1 and that obtained with Eq. (13), amount . 1, [se 
= Fic. 1. Body configuration. functi 


ing to 32 per cent atyw = 0.3 
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READERS’ 
x, and all other A and A constants given here are defined by 
Syvertsor 
(3) Cone-Cylinder-Frustum 
On an expanding conical frustum following a cone-cylinder 
iCy/da dCy/da)t7¢: (ARB? l 
1./Ko2) {(1 — e *2"/) ((1/2)Ke2 cot 07 + Koy + 1 
KL (6 
la - (1/3) (dCy/da A(2Re 1)? xX 
{ 
BY, + 2L, + (3/2) ((2N3/L,) + 1] cot 6; 4 
: Kol 9 - 
1./K l " atin Ys + (1/2) (KeX2 + 1) cot 6, 
2/K e Keks ((1/2)K» cot 0f + Ko Ly + Xs) + 2]} (7 
id 
( ( {R,? l 
€4( Cos ( tan 07!/Ke} [A 
2 i] (Ko + 2 tan 9, + 2L;/K,» tan 6;)e *? Ss 
re 
OP /Po)/OS 
K. = Q 
(P/F (P,/Po)| cos @ 
8 tan 67} tan 6,-(dCy/da (P,/P. FI A3{ 10 
e effective pitching moment, including the axial force con 
tribution, 1s 
1Cy/da - dCy da fy sec? 6, 
XY, tan 6, R y (dCx da tan 6, 11 
(4) Cone-Cylinder-Frustum-Booster 
On a evlindrical booster following a cone-cvlinder-frustum 
. Lp 
d( \ da B= \ | AK l é K3LB T A,Lp (12 
sldute = ~Aj/Ed 1M, + 0K 
¥, + (1/Ks) + Lele *2"9} A,Lpl|X5 + (1/2)La] (13 
ere 
{ SR (Ag / As) A l4 
1, = SRx (P./P Ng /Xd P./P. GA 15 
K; = [0(P/Po)/oS l P,/P 16 
(5) Cone-Frustum 
On a frustum following a cone, Eqs. (6) to (8) are applicable 


when the following changes are made 





1, replaced by A A» replaced by A, 
\ replaced by / Ci repl iced by C. 
8 tan 6,;|tan @--(dCy/da ANz/A,)A ‘ 
O(P/Po)/OS 
KK, = : . 18 
/ | P;/Po)| cos 6, 
(6) Cone-Frustum- Booster 
On a eylindrical booster following a cone-frustum, Eqs. (2 
nd (3) are applicable when the following changes are made 
1 replaced by A K replaced by A 
replaced by Lyx L. replaced by Xs 
ere 
Ag = 8Rp(du/dx)As 19 
WAP/Po)/AS} » 
K; = (20 
l P,/P. 

All of the 4, and A, constants are determined from relations 
which are defined in reference 1. The A, constants are related 
nly to the pressure discontinuities and pressure gradients at the 

tners, and are obtained directly from a form of Equation 9 of 


reference | since 


K = n/(X Xn 


Likewise, the 4 constants are related primarily to the loading, 


\, |see Equations (19), (B-22), and (B-23) of reference 1] and the 


function \, which is defined as 
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FORUM 


AX = 27p/(sin 2p 


where yu is the Mach angle It should be noted that the subscript 
tf refers to properties of a cone tangent to the frustum and travel 
ing at the same free-stream Mach Number 

Once these constants are determined, it is a relatively simple 
matter to obtain the incremental normal force and pitching mo 
ment of each body section from the equations presented herein, 
and the total characteristics are obtained by proper combination 


of the component properties 


REFERENCES 
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The Magnus Effect at High Reynolds Numbers 
Howard R. Kelly* and Ray W. Van Akent 
U.S. Naval Ordnance Test Station, ( 
Received June 19, 195¢ 


hina Lake 


A BRIEF DISCUSSION on negative Magnus forces has recently 
been published by Krahn,'! in which a physical mechanism 
for the production of the negative force is described, but no 
references are given to actual data It is the purpose of this 


note to review briefly important experimental results of the past 
and to present some current data that indicate progress toward 
Magnus fe 


a complete understanding of the ree on spinning 


evlinders (see Fig. 1 
The 


were probably those 


earliest quantitative force measurements evlinders 


of Laf i\ 


ir reproduced here sec 


His results for a smooth evlinder 


Fig. 1) as lift coefficient (¢ based on 





plan-form area, vs. the ratio of cylinder wall speed V to air 
speed [ rhese results do not correspond to a two-dimensional 
evlinder, since the evlinder used by Lafay had a fineness rati 
* Head, Aerodynamics Branch 
’ Aeronautical Engineer 
(POSITIVE) 
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Fic. 1. Magnus force on spinning cylinder and Magnus forces 
coefficients measured by Lafay 
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Fic. 2. Magnus force coefficients at various Reynolds Numbers 


of about 4, with no provision for eliminating end effects. How- 
ever, an indication of the effect described by Krahn is found in 
data available in 1910 

Results of further measurements made by Betz* and Thom‘ 
were similar to those of Lafay at the lowest Reynolds Number 
he used. We now believe that the critical Reynolds Number 
R. for the experiments of Betz and Thom was probably higher, so 
the ratio R/R, was low and no negative force would be expected 
The critical Reynolds Number is here defined to be the Reynolds 
Number at which the drag coefficient of a nonrotating cylinder 
decreases abruptly, due to the difference in flow separation with 
laminar and turbulent boundary layers. In the Betz experi- 
ments end plates were introduced, showing that the maximum 
lift coefficient at high spin could be increased to about 10, as 
compared to 3 or 4 for short cylinders without end plates. No 
one has yet measured a lift coefficient approaching the theo 
retical limit of 47 

Experiments with spinning cylinders have recently been com- 
pleted at the Case Institute of Technology in Cleveland. Al 
though the results have not been published as yet, the authors 
have learned that these experiments confirm the jumps pre 
dicted by Krahn and Lafay (see Fig. 1). A much closer approxi 
mation to infinite-cylinder conditions was obtained by measuring 
forces only on the center section of a 3-section cylinder. 

A program is now in progress at the Naval Ordnance Test 
Station to study the Magnus force on spinning cylinders over a 
wide range of Reynolds Numbers, both above and below the 
critical Reynolds Number. A sample of preliminary results is 
shown in Fig. 2. The tests are being conducted in the Convair 
subsonic wind tunnel in San Diego. The results shown are for 
smooth cylinders 30 in. in length and 4, 6, and 8 in. in diameter. 
In order to minimize the effect of finite length, disc-shaped 
end plates rotating with the cylinders were used. They are very 
much like those used in the experiments of Betz. Spin rates 
were varied from zero to 8000 r.p.m., and air speeds from 50 ft 
per sec. to 300 ft. per sec. allowed a variation of Reynolds Number 
from 1 X 10° up to 12 * 10°. The tunnel has a low turbulence 
level, as the measured critical Reynolds Number at zero spin 
proved to be about 3.5 & 10° 
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Only a qualitative agreement is found between the result 
and the analysis of Krahn. Even discounting the difference j 
critical Reynolds Number, the exact shape of the shaded regio; 
in Fig. 2 of Krahn’s article is not verified This is to be 
pected, since the transition and separation of boundary layer 
even on nonspinning cylinders is only partially understood, 
the presence of spin introduces several complicating factors 

The principal new result shown in Fig. 2 here is the contrag 
between results at very low and very high Reynolds Numbers a 
low spin rates. Apparently the laminar boundary layer is my 
less efficient in vorticity transport than the turbulent boundar 
layer. Physically, this is quite reasonable 

These tests are being continued in order to obtain results tha; 
are truly characteristic of infinite cylinders and to investigat 
the effect of critical Reynolds Number more thoroughly. Tw 
cylinder lengths, 30 and 60 in., are being used to evaluate and 
correct for end effects. Turbulence screens are to be inserted at 
various distances ahead of the model to obtain a much lower 
critical Reynolds Number. It is hoped that the ratio R/R. can 
be increased to such a point that a maximum lift coefficient 
(which usually occurs at about V/l’ = 4) can be obtained with 
turbulent boundary laver on both sides of the cylinder. It js 
possible that this maximum may be well over 10, even approach 


ing the ideal value of 4z 
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An Experimental Investigation of Flow over 
Blunt Nosed Cones at a Mach Number of 5.8* 


William T. O’Bryantt and Reginald M. Machellt 

Guggenheim Aeronautical Laboratory, California Institute of 
Technology, Pasadena, Calif 

July 2, 1956 


INTRODUCTION 


t 


measurements on a 40° half angle spherical nosed cone at 
Mach Number of 5.8. The purpose of the present investigation 


was to obtain more extensive information on the surface pressur 


i A RECENT EXPERIMENT Oliver? obtained static pressure 


distribution and the shock wave pattern on spherical nosed cones 
in hypersonic flow. Six models with cone half angles of 10°, 2 
and 40°, and bluntness ratios of 0.4, 0.8, and 1.064, were tested 
at a Mach Number of 5.8, the bluntness ratio being the ratio 


the radius of the spherical nose to the radius of the base of the 


) 


model. Static pressure measurements and schlieren observations 
were made on all six models at zero yaw and on two models at 
angles of yaw of 4° and 8 In analyzing the experimental re 
sults it was found useful to compare the actual pressure distribu 
tions with the predictions of the modified Newtonian approx! 

ce 


mation, C, = C cos? », where C, is the pressure 
x 


“max ma 

* This investigation was part of a broad hypersonic research program 
sponsored jointly by U.S. Army Ordnance and the Office of Ordnance Re 
search under Contract No. DA-04-495-Ord-19 The advice and guidance 
of Prof. Lester Lees in connection with this investigation are gratefully 
acknowledged 

+ Commander, U.S. Navy, under the sponsorship of the U.S Naval Post 
graduate School and the Bureau of Aeronautics, Department of the Navy 

¢ Lieutenant, U.S. Navy, under the sponsorship of the U.S Naval Post 
graduate School and the Bureau of Aeronautics, Department of the Navy 
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Se result efficient at the forward stagnation point, and 7 is the angle be 
ference tween th mal to the body surface and the direction of the 
ed region free-stream velocity. A complete report of this investigation is 
to's ae OB jven in GALCIT Hypersonic Wind Tunnel Memorandum No 
’ given in 
ry layer 39 
0d, ir 
re i RESULTS AND DISCUSSION 
| . 
contrast [he static pressure measurements on the two 40° half angle 
mbers at models and the three 20° half angle models at zero yaw are pre 
1s mu ented in Figs. 1 and 2. For both families of models the pressure 
yOundar distribution on the spherical nose agreed very closely with the 
modified Newtonian approximation, On the conical portions 
ults that of the models with half angles of 20° or more, the flow overex 
vestigat panded with respect to the Taylor-Maccoll cone pressure (as 
Vv. Tw tabulated by Kopal*) and then was recompressed to approxi 
late and mately this value The minimum pressure point moved down 
serted at stream as the cone half angle was decreased. On the conical 
h lower portion of 10° half angle model with a bluntness ratio of 0.8 the 
¢/ R. can werexpansion did not occur The similarity of the result 
efficient | | irs ) ilies ‘Is sl > ssur > 2 ° . - 
emicient within the first two families of models showed that the pressure Fic. 3 40° half angle spherical nosed cone at M 58. 
1ed wit! distribution was independent of the bluntness ratio, r/R, when Re/in. = 1.91 X 105 (schlieren photograph 
r. It is described in nondimensional coordinates. Further, the pressure 
yproach distribution was independent of the Reynolds Number over the 
range of test conditions 
- . r - — yy — 
A schlieren photograph of the flow over a 40° half angle spheri KEP o ie eines eine tout 
cal nosed cone at a Mach Number of 5.8 and a free-stream Rey | VilA ) AR : = 
Sciences, nolds Number per inch of 1.91 X 10° is shown in Fig. 3. Mid 8 oe OMA | at EXPERIMENT ” 
way out on the conical portion of the model the shock wave shows | Ly | RAY —_ 
s-Rendus ; a. ¢ er po — y . 
y se Cc re Ss . “te . Ss i “x , 
i slight reverse curvature his unusual characteristic was ob me ens oe Ps \\ 
. > ae ef pf “< 
tor. Zeit served only on the 40° models, and it was closely connected a = Fae . 
x ; : ‘ ‘ ‘ a hee: if } \A\\G 3 peo 
y, 1925 with the appreciable overexpansion and recompression which “omax oe aa g ~-a-o° A a om 
; Z 7 . wae pat 30 90 oskeaten 
ecurred on the conical portions of these models. From schlieren sais ——s | Yo ae 
Aer . ° . a eal es 
: photographs of the flow over the other models at Mach Number Zee | R0875 | oF aS ame 
120.35 fF hae 
2} i¢ at 
. { 
ROS J LOWER HALF ——+—— UPPER HALF fe— SCALE CHANGE 
wy oo cain cen © mmc? “30 20 10 08 06 04 02 0 02 04 06 08 iC : Fo 
‘ << S. ay R:0875 , = . : 
4 0 %)— : 
a \ ~ Prile2 1 os o4 238 - - 3 
an ——— EXPERIMENTAL Fic. 4. Surface pressure, vertical meridian plane, a GO”, 4°, 
| at NEWTONIAN 
over AN + ; — —— KOPAL & 
f 5.8* ' \, . 
° 
Comox \K + é ate 
TAA = i d = c ‘ 
4 a a 1 5.8 the separation distance of the bow shock wave from the 
ot Last “ . , 
LAST ORIFICE, ®2 | nose of the model was found to vary essentially linearly with the 
radius of the spherical nose 
Fig. 4 shows the effects of angle of vaw on the pressure dis 
tribution on a 40° half angle model, measured in the plane in 
r r : ; a 3 ¥. - which the model was vawed rhe upper half refers to the top 
ressure : half of the model when it was yawed above the free-stream direc 
e at : 7 P : : } 
aie ‘ } tion. If the angle of yaw is considered as a change in the ef 
— Fic. 1. Pressure distribution on spherical nosed 40° cones, jus? 5 ——— 
lea iy. - eo fective cone angle, then as the cone angle was reduced, the point 
Coeur of minimum pressure moved downstream on the conical portion 
Comte of the model, as in the case of zero yaw. On the spherical por 
yo tions of the yawed models the pressure distribution again agreed 
1 
testet t...7 CES RE ST SR neal Lhe the ead Ne J 
vy aT deli aes A teed very closely with the modified Newtonian approximation 
itio of nm 7 |, Smee MODEL EH Aelins i” i i ail seat | 
} * - |R: B75 3 «04 9 | Pressure drag coefficients obtained by integrating the unvawed 
f the : esl 86 banfay & 5 1064 19 eva : 
oo P iN f #3] 4) 4 08 3 pressure distributions agreed very closely with the predictions 
a \ tes” fea ft Srp : 4 08 238 : a : : 
ation ‘ sat 9? ____experiMENTAL of the modified Newtonian approximation, except for models 
lels at \ ia |. =a ‘th large cone rle 1 ll bl | he 
t 2c lee —-— KOPAL with large cone angies and sma juntness ratios, where the 
al re | } — | | drag coefficient approached the Taylor-Maccoll value for sharp 
tribu : i | nosed cones 
, Omron \ 
prox! x LAST ORIFICE, *5 
e ct "| \ aa aaa ams | REFERENCES 
¥ LAST ORIFICE, % ] 1 Machell, R. M., and O'Bryant, W. T., An Experimental Investigation 
ogram >| \ LAST ORIFICE, "3 Flow over Blunt Nosed Cones at a Mach Number of 5.8, GALCIT Hypersonic 
ce Re Xe g Wind Tunnel Memorandum No. 32, June 15, 195¢ 
ne ia ts tt 
idance ee 2? Oliver, R. E., An Experimental Inve tion of Flow about Simple Blunt 
tefully “— SCALE CHANGE Bodies at a Nominal Mach Number of 5.8, Journal of the Aeronautical 
a os ._ _ on 30 40. 50 Sciences, Vol. 23, No. 2, pp. 177-179, February, 1956 See also GALCIT 
| Post 3. Hypersonic Wind Tunnel Memorandum No. 26, June 1, 1955 
ivy . § Staff of the Computing Section, Center of Analysis, Mass. Inst. of Tech 
} ) 
Post IG. 2. Pressure distribution on spherical nosed 20° cones, under the direction of Zdenek Kopal, Tables of Supersonic Flow Around 


Vy Mo = 5.8 Cones, Technical Report No. 1, Cambridge, 1947 
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Some Comments on Surface Contamination 
in Hydraulic Analogy Research 


R.A. A. Bryant 
Lecturer in Mechanical Engineering, the New South Wales 
University of Technology, Sydney, Australia 


July 9, 1956 


i A RECENT NOTE Bomelburg! has stated that it is by no means 


an easy matter to carry out successful experiments using the 


hydraulic analogy; it was also stated that even the accumulation 
of dust on the free surface affects the reliability of the analogy 
more or less 

From experience in operating a similar shallow water tank the 
present writer agrees entirely with the first statement and has 
found that both dust and contamination of the surface of the 
model are sufficient in themselves to make experimental results 
quite worthless. However, a somewhat different approach to 
that of BOmelburg was made in overcoming the dust problem 
It was felt that water should be used as the test fluid since the 
use of kerosene means that one must either decrease the depth 
of the liquid in order to obtain surface wave propagation with 
zero dispersion (for all wave lengths greater than one inch) or 
operate with finite dispersion at about the same depth as for 
water. A reduction in depth is not desirable as it further com 
plicates an already difficult measurement problem; likewise the 
latter alternative is undesirable as it produces a partial break 
down of the analogy 

In early experiments it was found that dust on the surface 
had a twofold effect. It not only varied appreciably the surface 
tension of the water but also adhered to the surface of the model 
(brass) and produced nonuniform wetting 


1 


Consequently, a technique was developed to produce an un 
contaminated water surface; it involved careful cleaning of the 
tank prior to adding the water, after which the surface was swept 
with a perspex runner (see Fig. 1). The underside of this runner 
just contacts the free surface so that on traversing it along the 
tank the dust is carried before it and finally accumulated at the 
end for removal with a small pump. That such sweeping will 
cope with dust concentrations considerably in excess of those 
normally experienced is shown by Fig. 2 

Production of an uncontaminated water surface improved 
matters but nonuniform wetting of the model still occurred 
The inception of this was traced to chemical action taking place 
at the interface of the model and the water. During continued 
operation of the tank more and more of the model was affected 
in this way until the surface characteristics had become different 
from those existing originally and, according to the depth of 
submergence, a band of surface existed having variable surface 
activity along the whole length of the model. On treating the 
surface of the model with silicone fluid the chemical action was 
stopped. Uniform nonwetting conditions were then attained 
which remained until dust, precipitated from the air, had again 
contaminated the water surface. However, because of the high 
viscosity of the silicone, only very minute amounts of dust had 
to be present for a considerable accumulation to occur on the 
model. The overall! effect was worse than for the untreated 
model, and, although the water surface could be quickly cleaned 
the cleaning and recoating of the model were both time-consuming 
and awkward 

It was finally found that the best coating for the model was 


a clear resin-base lacquer (clear ‘‘Dulux”’ paint) which is non 


wetting. The lacquer was applied by means of a steel roller in a 
number of very thin coats, until a thickness of about 0.01 in 
had been built up. Such thickness was found necessary for the 
grid scribed on the side of the model to be completely sub 
merged and the surface to be without discontinuities 


It has been found that, by treating the model in this way and 
sweeping the water surface, the tank may be operated for periods 
of about one hour before the concentration of dust becomes 


significant. A resweeping of the surface may be carried out in 
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Fic. 1 Perspex runner. Rubber strips at either end form a s 
with the vertical sides of the tank 











Fic. 2 
taken normal to water surface from above 
below glass bottom of tank ; 

Top. Direction of sweeping from left to right The few 
particles to the left of the runner were grains of sand resting 
the glass plate, which were accidentally included with the dust 
The dust was produced artificially in order that clear photograp® 
could be obtained ; 

Bottom. Result obtained by reducing contaminated area fro! 
176 in. by 39 in. to 1'/2 in. by 39 in 
here produced between two runners but is normally isolated at 
one end of the tank 


Effect of sweeping the water surface. Photographs 
Illumination [rot 
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minutes, while the surface of the model is inert and 


less than five 
requires n ttention 
REFERENCI 
Bome!lbur H. J.. Comment on Wedge Pressure Coe ficients in Transont 
by H ‘ {nalogy. Journal of the Aeronautical Sciences, Vol. 22 
731. October, 1955 


Errata—‘‘Minimum Weight Design of a Simply 
Supported Transversely Stiffened Plate 
Loaded in Shear’’! 





Malcolm F. Symonds 
uctures Engineer, Hiller Helicopters, Palo Alto, Calif 
August 2, 1956 
l In the Symbols, the ¢ in “¢ = average plate-stiffener 


tickness”’ should be 7 


2) Fig. 1 and Fig. 3 have been incorrectly transposed 


3) In Eq. (12), “*(g/h ’ should be ‘‘(q/h)/7 


REFERENCE 
Malcolm F., Journal of the Aeronautical Sciences, Vol. 23, 


No. 7, pp. 685-693, July, 1956 


Symonds 


Propulsive Ducts at the Tips of Rotor Blades 


P.R. Payne* 
Bristol Aero-Engines Ltd., Filton, Bristol, England 
uly 9, 1956 


A*°™ presented by Foa, Gail, and Goodman! in the April, 
Readers’ dealt with the problem of the 


optimum cross-sectional 


1956, Forum 


structurally shape for a rotor tip 
mounted ram-jet tailpipe 

The use of a flexible tailskin was first suggested by the present 
writer in reference 2, and is the subject of a Patent specification.’ 
The idea was first used in practice by the writer early in 1954, at 
\uster Aircraft Ltd., (England) and later at Bristol Aero-Engines 
Ltd 

\n analysis similar to that of Foa, Gail, and Goodman was 
made in 1953, but during the more detailed work of designing 
the first ramjet it was realized that the two-dimensional approach 
was an oversimplification. Although the thin skin is unable to 
carry bending moments in the plane of the paper (see Fig. 1 of 
reference 1) the aggregate loading on the other elements results 
in shear stresses which are not considered in reference 1 In 
practice it is essential to relieve this shear by the provision of a 
“beam” inboard of the main body of the ram-jet, to which the 
centrifugal loads of all cross-sectional elements are transmitted. 
If this beam is not inside the tailskin shell, the point of maximum 
stress is at the joint Q in Fig. 1 of reference 1. In many practical 
applications it is not possible to make this joint strong enough if 
the tailskin is of uniform thickness, so that a preferred solution 
is to mount the beam inside the tailskin, such that the tailskin 
“hangs” from the beam under the influence of centripetal acceler- 
ition. If the beam is deeper than that required to fit the ideal 
shape of reference 1 the hoop tension in the tailskin will force it 
to fit tightly to the beam. Much of the tension is then diffused 
to the beam by friction, so that the joint Q is relatively un- 
Stressed. 

This method has been used in practice up to tip Mach Numbers 
of 0.8, the hollow beam being cooled by air tapped from the ram- 


ia : p : 
jet diffuser, with results that have been eminently satisfactory. 


* Now, with A. V. Roe Ltd., Malton, Ontario, Canada 
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Moreover, the finite depth of the beam at the exit plane of the 
ram-jet reduces the tailskin stresses in the most critical region 
Foa, Gail, and Goodman do not appear to have realized that 
there is no internal overpressure in the exit plane of a subsonic 
ram-jet, a fact which completely invalidates the two-dimensional 
approach 

In conclusion, it might be as well to point out that the thrust 
of a ram-jet is governed by the exit area, and not the maximum 
area as implied in reference 1. Moreover, minimum periphery 
is not an over-riding requirement, since, although it reduces flam« 
chilling and skin friction, it is easier to achieve satisfactory flam« 
spreading in a more elliptical section, because of the limitations 
imposed on baffle shape by the combination of high temperature 
and high centripetal acceleration. It should also be realized that 
a lower limit to skin thickness is set by the pressure pulsations 
which occur in even the smoothest burner, so that in practice it 
an equivalent ellipse ratio of less 


1S rarely possible to achieve 


than 2:1 
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On the Flexure of a Stepped Cantilever Beam 


E. Edryd Jones 

Department of Mathematics, The University of Nottingham, England 

July 13, 1956 

B* AN EXTENSION of a method due to Brown,! and use of the 
Laplace transformation, the flexure of a single-step canti 

lever beam with uniform load has been investigated by Carslaw 

and Jaegar.2 However, the resulting expression for the deflec 

tion of the beam becomes unwieldly 

continuities of the flexural rigidity of the beam is greater than 

one. The flexure of a nonuniform beam, with prescribed end 


when the number of dis 


conditions and loads, has been investigated by the author, 
and in this note the methods of references 2 and 3 are combined 
to determine the deflection of a stepped cantilever beam in a 
form which is readily adapted to numerical interpretation 

In the theory presented here it is assumed that a beam AB 
is clamped horizontally at the end A, (0, 0), and that the other 


end point B at x = L is free, with the y-axis directed vertically 
downward. The beam has a flexural rigidity B(x) at the 
point (x, y), and the bending moment there is M(x) = B(x) X 
d*y/dx*. The cantilever beam is unique in that the conditions 


it (0, 0) concern y and its first derivative, while those at x = 
L concern the bending moment and its first derivative, and these 
four conditions can be introduced into the solution of the equa 
tion of flexure by considering the bending moment at station x 
due to the loading from that point to the free end of the beam, 


cf. reference 2 


It is assumed that if 7(p) is the Laplace transform of y(x), then 
ec =x 
Wp) = | e~”* y(x)dx 
J 
The bending moment equation is d?y/dx? = B(x)M(x), and its 
Laplace transform is 
p? ip) = BM (1 


since y = dy/dx = 0 when x = 0 
considering its inverse by use of the convolution integral (see 


Dividing Eq. (1) by p?, and 


reference 2, p. 80), then 


W(x) = [ (x 7) B(x 
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/ 
If the prescribed bending moment at B is Wo, and the shear force The Effects of Prandtl Number on High-Speeg | repla 
5 ° | 
is M,, then Viscous Flows over a Flat Plate* 
M(x) = Myo + ML — x) + m(x (3) 
Y.H. Kuo 
where m(x) is the moment due to the loading on the beam be Professor, Graduate School of Engineering 7 
tween station x and B. A substitution of Eq. (3) into Eq. (2 Cornell University, Ithaca, N.Y | “ 
gives the deflection of the beam in terms of the given end condi July 16, 1956 sti 
tions mp! 
| 
: HE FLOW of a viscous perfect gas over a semi-infinite 
ISX AMPLE , ; : 
plate at high flight Mach Number J/ was considered jn 
The example considered here concerns the flexure of a stepped previous article! under the assumptions that (1) the plat 
cantilever beam under its own weight. The problem is compli- thermally insulated, (2) the Prandtl Number Pr is equal 1 ae 
cated by the fact that both A(x) and the loading function w(x) unity, and (3) the viscosity w varies linearly with temperature 7 a 
are discontinuous at points of the beam If there are discon With the choice of a special coordinate system, the problem wa It t 
. ens : , ° . ‘ i ‘ : Macl 
tinuities at the s points x = x,, then the following functions can treated by a method of perturbation, with the case of unifor; Mac 
be formed: yressure as the zeroth approximation. To the first order, t 
I PI : 
i “te method vields a velocity field and pressure which approac! 
B(x) = z a w(x) = a Wy ; , ; hs PI “ae 
ead aes large distance from the plate the simple wave solution. Being —" 
of the simple wave type in the inviscid field, these solutions 
> > o are Cc ‘ _ ») ‘ . = 
where b,, w,, are constants, and n,m = 0,1,2,..., 8, with 0) actually correct to second order there 
re H a 7 -) iia Lor i 7 2%) oc ini on . . . : 
With this form of B(x), the integral in Eq. (2) containing The results of this analysis reveal an important parameter for th 
and \/) can be written as (defined further below), besides the Reynolds Number R t can 
b 6 Fa fy WL ’ (= U,L/ve, Uo, vo and L being, respectively, the velocit SONIC | 
’ (x v) + ML — v)j di : ; Maes os : 
r bn : cea sé and viscosity of the gas at infinity) Pogether they appear it lows 
il all flow quantities. Unlike the Reynolds Number, which char 
which has the value acterizes the viscous effects, v depends primarily on the ten 
1 ke perature distribution and, therefore, is a function of Mach Num for the 
} . , \2 {2 } 3 — Ir -\t og. ogc 
B secx One Xn)? (SM M(3L sade = (4 ber and represents the compressibility effects of the gases. It 
; : eee ; is the purpose of this note to investigate the functional depend 
Phe function m(x) is given in general by ence of v on both Mach and Prandt! Numbers for both the case — for hee 
E of an insulated plate and that of heat transfer—i.e., cold wall it high 
mx) = (u — x)w(u)du As shown in reference 1, the boundary layer deflects the mait reduce 
sii flow at a given point to the first approximation, by an amount Numbe 
and on integrating by parts which is proportional to z § insulat 
Pris a 
1 a. % : . ° 
mx) = 5 (L xPw(Ll)—- 5 £4 (%m — x)? w, (5 vw = 0.860 tT’ dt l ence, 
“« = Im>zT 0 
; ' In t 
since w(“) changes by w,,, when u passes through u = x, Sub where 7) is the temperature ratio corresponding to constant | bow 
stituting from Eq. (5) into the integral of Eq. (2), which contains pressure in the main stream and ¢, the similarity variable. Th ransfe 
m(x), then the contribution to the deflection is prime under the integral sign denotes total differentiation with | eee? 
os respect to ¢. In the case of constant Prandtl Number an Leware 
: w( L) | (x — v) (L — v)2B(v)dv — constant wall temperature, the temperature 7) was found tobe Fy Joel] 
— T, = 1+ [(3 1)/2] M? @0(¢, Po 2 
erm = 9 
l . ; 2 
= ym Wm (x — v) (Xm — v)*B(v)di for the insulated plate and 
= Im Zz 0 . 
een F . > . . \u 
sa To = 1 + [(y — 1)/2] MP? Ow(E, Pr) 4 
I b» W, (x — v) (Xm — v)*B(v)di 41 + [(y — 1)/2] MeV" Pr — T,,} On0(¢, P 
2 tm>z 
Ji) where j 
: : ; y : : for heat transfer. Here 7, and y are, respectively, the wall tem parame 
On substituting the step function form of B(v), and evaluating perature ratio and the ratio of specific heats. The functions ing edg 


simple integrals, the contribution reduces to 6:9 and 6 are defined in terms of the Blasius function fi(f 
l “— : For both cases, vp as an explicit function of y, 17, T,, and Pr cat 
D4 w(L) > Ds bi(L — x)t + (L — x)? (4x — L — 3x,)} - I 


2 In<i be found by evaluating two integrals 


ss cane fi | 
WmDn (Xm — Xn)? (4¥ — Xm — 3xXn) — = vf f."\PT pm \2—Pr 7. . | 
24 rm<cx tn<Im es . i= x Uo ) A ( di di 


l , w\ Pr Jie 
37 p bd WmOn \(Xm x)t + (xm — Xn)? (4x — Xm — 3xn)f pe a f . — ' 
24 rm>z rn<7 and I, = fi ¢ dt ) — and, fin 

¥ 


(6) 


The complete deflection is given by Eqs. (4) and (6), and where ap denotes the initial value of f ( 
can be used, for example, to determine the effect of the presence Numerical integration shows that J; is 1.54 and 2.04 for fy | 
of notches in the beam. If the loading is uniform, so that equal to 0.72 and 0.50, respectively. From this, together wit | ¥ assu 
wo = w(L) = W,w, = Oforn > 1, then only the first series term the value at Pr = 1, it is proposed to approximate /; by the | where ; 
of Eq. (6) remains formula T* In 
The analysis may be extended to cover point loads and moments lL, = 1.192 (Pr)~°7 t = 
A i . ali a | 
and distributed moments 
This agrees with the exact values at Pr = 1 and 0.72 and gives4 |, 
REFERENCES value of 2.00 at Pr = 0.5, with an error of only 2 per cent. It i 
> cace of > > “OXI i ly to 
Brown, C. I The Treatment of Discontinuities in Beam Deflection the case of J», the usual method of approximation, name It is 
Problems, Quart. Appl. Math., Vol. 1, p. 349, 1944 . , &reatly 
? Carslaw, H. S., and Jaegar, J. C., Operational Methods in Applied Mathe * This work was carried out as part of a research program at the Graduat tically y 
matics, 2nd Ed., p. 327; Oxford University Press, 1948 School of Aeronautical Engineering and was supported by the Office of Nava lecr 
7 ; lecrease 
Jones, E. E., The Flexure of a Non-Uniform Beam, Pac. J. Math., Vol Research Reproduction in whole or in part is permitted for any purpoe : ' 
luction « 


5, Supp. 1, p. 799, 1955 of the United States Government 
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replace y) in the integrand by exp | —aet#/12], leads to 
(- V3 3r[l1 + (2/3 4.93 
= = (7 
3 (aoPr)? (Pr)? 

On the ind, at Pr = V.72, numerical integration yields an 
iswer of 6.55, while the approximate Eq. (7) gives 6.14. To 
iprove this approximation, Eq. (7) is modified to read 

Ts, = 5.182/(Pr)? (8 

This formula again agrees with the exact value at Pr = 1 but 

fives a Value 6.45 at Pr = 0.72, which is in error by about 1.6 per 

g 
ent 
If the approximate Eqs. (6) and (8) are accepted, at very high 


Mach Number zp can be approximated by 


0.596 (> 1)Pr Vr? (9 
for the insulated wall; and 
0.596 Pr 0.430Pr'6] (5 — 1)M (10 


for the case of coo From these 


it can be concluded that for Pr different from unity the hyper 


ng with low wall temperature 


sonic parameter X (see reference 3) has to be modified in the fol 


lowing manners It is 
Pr V3Re 2 (11 


for the insulated plate; and 


PA-.2% — (2.721 Pr M? Re~‘"/2 (12 


for heat transfer. It is further noted that, in the case of cooling 
it high Mach Number, the compressibility effects will be greatly 
reduced, in fact, the induced pressure at the same Reynolds 
Number is almost 70 per cent less than that if the plate were 
insulated. In the case of the insulated plate, the effect of reduced 
Pris a slight reduction of the deflection of the streamlines and, 
hence, an attendant reduction of the induced pressure 

In the problem of high-speed flight, it is of importance to 
heat 


know, in addition to the local friction, the rate of local 


transfer at a specified constant wall temperature. For arbitrary 
Prandt! Number, a calculation similar to that in reference 1 can 


be carried out. This gives a local skin-friction coefficient c; and 


ilocal heat-transfer number .Vu, as follows: 


= 2V « Rew P/é a Vv < Re (vay BV &) + (13 
Nu = P) \ cRe V p & |a@ Vv c/Re (1 a) BY é + 

(14 

where p and & are, respectively, the local pressure ratio and a 


parameter, both being functions of the distance x from the lead 


ing edge;' a, a;, and 6 are defined by 


a = (ao/2 4 M2 4 (y + 1)/2] (4/82) — 2 — (w/a { 
sy + 1)M4 24M? ' 
a =a 132 
Pe ar (15 


j 


ind, finally, the Chapman-Rubesin constant C, in the present 


Case, is Chosen as 
C = [0.45 + 0.557, + 0.09(7 - 


1)M24/ Pri®-! (16 


by assuming that the viscosity-temperature law at the point 
where it coincides with the linear relation is proportional to 
I. In the case of large M, a and & can be approximated by 


) 99) 


0.596Pr*" 4 
0.480Pr''6} (4 1)M?}> (17 
$] 2 \ 


It is interesting to note that both constants a and % are 
gteatly affected by cooling but the product a;v» remains prac- 
ucally unchanged by heat transfer. Thus, insofar as cooling can 
lecrease the skin friction, it is effected mainly through the re 
luction of pressure 
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Comments on ‘‘The Ultimate Strength of 
Multiweb Box Beams in Pure Bending’”’ 


Richard A. Pride 


Aeronautical Research Scientist 
NACA, Langley Field, Va 
July 27, 1956 


Structures Research Divisior 


oo NEEDHAM! PRESENTED an empirical method for 
predicting the ultimate bending strength of multiweb box 
beams having formed sheet-metal webs. On the basis of the 
agreement with the results of 49 tests,” a reader justifiably could 
conclude that the method adequately accounts for the factors 
which determine the strength of this construction. This reader, 
however, believes that the applicability of the method should be 
more carefully qualified inasmuch as it is based upon a test series 
in which certain riveting parameters were assigned special values 
These parameters are not included in the method and can have a 
marked influence on multiweb beam strength as will be shown in 
this note 

The 49 beams in the experimental program? had very similarly 
skin. The 


along the web attachment flanges had a pitch equal to three rivet 


designed riveted joints between webs and rivets 
diameters, and their centerline was offset from the web plane a 
minimum possible distance consistent with a 4 fy bend radius on 
the flange. This particular rivet pattern was selected because in 
creases in rivet pitch are known to decrease the crippling strength 
of sheet-stringer panels,’ and a preliminary multiweb-beam test 
program‘ had shown that increasing rivet offset caused a decrease 
in strength. The riveting done by air-frame fabricators is sub 
ject to considerable variation from the specifications of this test 
series. There were three variations in rivet pitch and offset 
associated with the three skin thicknesses in the test beams 

The failure criterion developed by Needham is based upon the 
assumption that beam failure is associated with the crippling 
stress for the attachment flange, and he chose the width-thickness 
ratio of this flange as being the key parameter. His only re 
quirement as regards riveting is that the pitch shall be small 
enough to prevent interrivet buckling and that the rivet-line 
offset shall be a minimum consistent with good fabrication prac 
tice. In the test beams there were three changes in flange width, 
but it is important to note that the ratio of rivet-line offset to 
flange width was fairly constant 

A good check on the method is provided by tests of beams which 
have varying ratios of rivet offset to flange width. Consider the 
four beams reported in the preliminary multiweb-beam test pro- 
gram.* These beams were fabricated of 7075-T6 aluminum alloy 
and were nominally identical, except for changes in the bend ra- 
dius of the web attachment flange which led to variations in the 
flange width and rivet offset. Measured dimensions and failing 
moments are given in Table 1 (nominal dimensions were given in 
reference 4). The beam with a flange width of 1.030 in. was 
fabricated using angle extrusions for the webs rather than formed 
sheet metal; however, this should make no basic difference in the 
calculation procedure. It can be seen that Needham’s method 
calculates the failing moment of the beam with a flange width of 
1.098 in. very well, but it is too low for the beams with smaller 
flange widths and too high for the beam with the larger flange 
variation in beam 


width. This comparison shows that the 


strength is greater than predicted on the basis of the variation in 
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TABLE 1 


Experimental Data for 7075-T6 Aluminum-Alloy Formed 
Channel-Web Multiweb Beams with Variations in Bend Radius 


Allbeams: ov, 2.2 ksi except webs of beam with 1.030 flange 
width, o-, = 78.0 ksi, E = 10.5(10)% ksi, rivet diameter = 
5/16 in., rivet pitch = 15/16in 
Flange width, } 1.030 0.968 1.098 1.228 
Rivet-line offset 0.44 50 360.65 ~=—0..74 
Bend radius 0) 2/16 1/16 6/16 
ts 0).239 249 0.247 0.249 
tu 0.060 O64 063 0.064 
Cell width, bs 7.50 . 50 7.47 7.48 
Beam width, B 24.14 OO 24.01 24.00 
Beam depth, h 5.68 67 5.69 5.72 
Experimental failing moment, 
in.-kips 
Calculated failing moment, in.- 
kips (Needham's method ) , »2e 1,145 
Deviation from experiment, % 3: 2 +10 


.400 ,220. 1,045 


flange width and suggests that the much larger variation in rivet- 
line offset in these test beams is responsible for the trend in beam 
strength 

A recent investigation’ reported the results of 16 additional 
multiweb-beam tests in which 7075-T6 aluminum-alloy cover 
skins were stabllized by channel-type webs fabricated with ex- 
truded angles. These beams had rivet-line offset distances much 
smaller than any beams previously discussed and also had smaller 
ratios of rivet offset to flange width. Needham’s method under- 
estimates the strength of these beams by an average of 18 per cent. 

In order to indicate the influence that variations in rivet pitch 
may have on failure strength of multiweb beams, results for three 
tests not previously published are given in Table 2. These three 
beams were fabricated to be nominally identical except for varia- 
tion in rivet pitch. Changing the pitch from 4 rivet diameters to 
12 rivet diameters resulted in an 18 per cent decrease in experi- 
mental failing moment for this particular beam cross section. 
It should be noted that the 12 diameter pitch (2'/; inches) is more 
than adequate to prevent interrivet buckling—the only restriction 
placed on rivet pitch by Needham. Needham’s predictions are 
considerably higher than the experimental failing moments of 
these three beams and, of course, indicate no variation due to 
change in rivet pitch. The slight variation shown in the three 
predicted moments result from variations in actual beam dimen- 
sions as indicated in Table 2. The high overall predictions can 
be traced to the relatively smaller attachment flange width-thick- 
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TABLE 2 


Experimental Data for 7075-T6 Aluminum-Alloy Formeg 
Channel-Web Multiweb Beams with Variations in Rivet Pitgh 


All beams: o,, = 72.2 ksi, E = 10.5 (10) ksi, rivet diameter = 


”] 


3/16 in., attachment flange bend radius = 4ty 


Flange width, b 0.78 0.78 0.78 
Rivet-line offset 0.50 0.50 0.50 
Rivet pitch 9 1 
ls 0.250 0.250 0.247 
tn 0.0658 0.0670 0. 0662 
Cell width, bs 7.50 7.49 7.50 
Beam width, B 23 .33 23 .32 23.30 
Beam depth, h 5.29 ae 5.26 
Experimental failing moment, 

in.-kips 1,065 875. 
Calculated failing moment, 

in.-kips (Needham’s method 1,259 oe 1 , 229. 
Deviation from experiment, % +18 


" , ‘ 
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ness ratios in these beams than in comparable beams of references 
2 and 4 

From the above comparisons it can be seen that small varia. 
tions in the offset of the rivet centerline from the web plane and 
variations in rivet pitch produce significant changes in beam ulti- 
mate strength. Such variations are inevitable in beams of dif- 
ferent design and are well within the limitations that Needham 
placed on the method. It can be concluded, therefore, that, as 
an empirical method based on tests of beams with a particular 
riveted attachment flange design, the method will give uncertain 
results when applied to beams of slightly different detail design. 
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